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FOREWORD 


This report presents results of the expansion and improvement of the 
FORMA system for response and load analysis. The acronym FORMA stands 
for FORTTIaN Matrix Analysis. The study, performed from 16 May 1975 
through 17 May 1976 was conducted by the Analytical Mechanics Department, 
Martin Marietta Corporation, Denver Division, under the contract NAS8- 
31376. The program was adniinistered by the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, Huntsville, 
Alabama under the direction of Dr. John R. Admire, Structural Dynamics 
Division, Systems Dynamics Laboratory. 

This report is published in seven volumes: 

Volume I - Programming Manual, 

Volume IIA - Listings, Dense FORMA Subroutines, 

Volume IIB - Listings, Sparse FORMA Subroutines, 

Volume lie - Listings, Finite Element FORMA Subroutines, 

Volume IIIA - Explanations, Dense FORMA Subroutines, 

Volume IIIB - Explanations, Sparse FORMA Subroutines, and 
Volume me - Explanations, Finite Element FORMA Subroutines. 
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ABSTRACT 


This report presents techniques for the solution of structural 
dynamic systems on an electronic digital computer using FORMA (FORT RAN 
tfatrix ^alysis) • 

FORMA is a library of subroutines coded in FORTRAN IV for the effi- 
cient solution of structural dynamics problems. These subroutines are 
ir the form of building blocks that can be put together to solve a large 
variety of structural dynamics problems. The obvious advantage of the 
building block approach is that programming and checkout time are limi- 
ted to that required for putting the blocks together in the proper order. 

The FORMA method has advantageous features such as: 

1. subroutines in the library have been used extensively for many 
years and as a result are well checked out and debugged; 

2r method will work on any computer with a FORTRAN IV compilei:; 

3. incorporation of new subroutines is no problem; 

4« basic FORTRAN statements may be used to give extreme flexi- 
bility in writing a program. 

Two programming techniques are used in FORMA: dense and sparse. 
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LIST OF SYMBOLS 


[ ] 
{ } 

{ f 

T 


matrix 

column matrix ^ 

> vec tor 
row matrix ^ 

transpose (when symbol is a superscript) 


JnP^n the row size of matrix 

n designates the column size of matrix 

a^^j a designates an element of matrix [a] 
i designates the ith row of matrix [a] 
j designates the i th column of matrix [a] 



I-l 


I. INTRODUCTION 


This volume presents an explanation of the function of each dense 
subroutine in the FORMA library. Example problems are given in some 
cases to clarify the operations performed by a subroutine. 
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II. SUBROUTINE EXPLANATIONS 

The subroutines are given in alphabetical order with numbers 
coming before letters. 



AABB 


Subroutine AABB calculates the summation of two matrices, 
each matrix multiplied by a scalar. In matrix notation, 

^^^NRxNC “ ‘**-^^NRxNC ®^®^NRxNC 


where 


aij = u a^^ + 3 b^j 


/i = 1, m\ 

U = 1. Nc/ 


NR is the number of rows of each matrix, and 

NC is the number of columns of each matrix. 

The number of rows of [A] and [B] must be equal, and the number 

of columns of [A] and [B] must be equal. 

Theorem : Matrix summation is commutative. 

That is, [C] + [D] = [D] + [Cj. 

Theorem ; Matrix summation is associative. 

That is, [C] + ([D] + [E]) * ([C] + [D]) + [E]. 

EXAMPLE 


Consider input of u = 5., 3=2., 

7. 2. -3. 

4. 5. 6. 

The reader can easily verify the output to be 


[A] 


2x3 


, and [B] 


2x3 


7. -8. 9. 

10 . 11 . 12 . 


tz] 


2x3 



’ 7 . 

2. 

-3.' 


' 7. 

-8. 

9.' 

- 5. 




- 2. 





_4. 

5. 

6.^ 


10. 

11. 

12. _ 


21. 26. -33. 

0. 3. 6. 



AIaHM--!/ 


Subroutine ALODl takes (on option) distributed and concentrate, 
lateral forces on a beam and replaces them with representative 
concentrated forces at selected points on the beam. 

The x-stations of the selected points (panel points) are given 
in {PP}. These x-stations must be in increasing order. 

The distributed lateral force, f(x), is assumed to be piece- 
wise linear and is represented by straight line segments as shown 
in Figure 1. 


f(x) 



Figure 1 Distributed La^'eral Force 

The x-stations of the end points for the line segments giving the 
distributed lateral force are independent of the panel point 
x-stations. However, the distributed lateral force must be within 
the panel point limits. The line segments representing the dis- 
tributed lateral force may or may not be joined and may overlap. 
The distributeu lateral force is defined in [DIST]. Each row cf 
[DIST] represents one nonvertical line segment. The form of each 
row of [DIST] is [x] x^ f 2 ] where xj, f^ give the first end 
point, and X 2 , ^2 8^ve the second end point of a line segment. 

The concentrated lateral forces are defined in [CONC]. Each 
row of [CONC] contains one concentrated lateral force, F^, and 

Its X location in the fonnjx^ fJ • A concentrated lateral force 

may be outside the panel point limits. 

The calculated representative concentrated forces at the 
selected panel points are placed in {Z}, The total lateral force 
and center of pressure are also calculated and printed. 



ALODl—2/7 


DESCRIPTION OF TECHNIQUE 

The replacement of distributed and concentrated lateral forces 
by representati\e concentrated forces at selected panel points is 
obtained using a virtual work approach. The virtual work done by 
the lateral forces on a beam is defined by 


6W 



f(x) 6y(x) dx + 




( 1 ) 


where 

f(x) is the distributed lateral force, 

F is a concentrated lateral force, 
c 

6y(x) is the lateral virtual displacement measured from 
the body x—axis , 

X is the undeformed longitudinal axis of the beam, 

Xg is the starting x-station of the beam, and 

x^ is the ending x-station of the beam. 

The finite summation is over the number of concentrated lateral 

forces, F . 

c 

Most techniques for. describing the lateral displacement 
(virtual or real) assume a function, y(x), over the entire length 
of the beam. However, considerable skill is required in choosing 
this function. The technique used here is to represent the lat- 
eral displacement between consecutive panel points by some simple 
function of the panel point displacements. A linear* displacement 
function will be assumed here. This is illustrated in Figure 2 
where Xj^ and consecutive panel point stations. The 

region between panel points k and k+1 is referred to as bay k. 



Figure 2 Linear Displacement Function 
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The displacement between panel points k ;md U+l is given by 
y(x) = + (^ - \) (Vk+i - yk)/(\+i - \)- 

Similarly 

6y(x) = + jx - x^j - x^) • 

The virtual work of distributed and concentrated forces will 
be considered separately. The distributed lateral force, f (x) , is 
considered first. The geometry for a line segment of distrib- 
uted force is shown in Figure 3. 


f(x) 



XI X2 

Figure 3 Line Segment Geometry 

li vi equation for a straight line segment as showr. In Figure 3 is 

f(x) = fi + (x - xi) (f 2 - fi)/(x 2 - xi). (3) 

Substituting Equations (2) and (3) into (1) gives the virtual 
work of the distributed lateral force represented by one line 
segment, i, in bay k as 
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The subscripts p aud q have been introduced to handle the possi- 
bility of a line segment extending past the bay limits. Thus* 

X is the greater of xi or x, and x is the lesser of x- or x, . , 

p * k q ' k +1 

Similarly, f is either fi or f, , and f is either £9 or f , . , • 

p ^ k q ^ k+l 

The integration is continued for the line segment »n adjacent 

bays, if necessary, until the entire lin^ segment has been used. 

Performing the integration of Equation (4) yields 



where 

^ ■ s ('r (“Ic+I - -p) - ‘-p} ^ (Vi - «p) - • 

- % ('p (% - \) + ‘•p} ^ («p - M * '‘■p})/“tc • 

L = X - X , and 

P q P 


•■k ■ - “k- 


By definition, the rroefficient of virtual displacement is the 
force. Thus, and ^re the representative concentrated 

forces at panel points k and k+l replacing the distributed lateral 
force represented by one line segment i in bay k. 


The concentrated lateral forces are considered last. For a 
force located in bay k, the displacement at the force location. 


oy 


is given in terms of adjacent panel point displacements 


(M' 

^y(*k+l)* from Equations (2) and (1) the virtual 

work of one concentrated force, F , in bay k is given as 


<SW 


c,k 




(5) 

(5a) 

C5b) 

(5c) 

(5d) 


( 6 ) 
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where 

\ (\+l “ 

and 

Vl ■ ("c ' ■ M' 

The virtual displacement coefficients, and the repre- 

sentative concentrated forces at panel points k and k+1 to replace 
the concentrated lateral force, F^, in nay k. 

The representative concentrated panel point forces for the 
entire beam are finally obtained by evaluating Equations (5) for 
each line segment of distributed force and Equations (6) for each 
concentrated force. Like terms of z are then summed. 


It is interesting to note that the virtual work approach with 
assumed linear displacement between consecutive panel points used 
here, gives identical results to the more common static "beaming- 
out” approach. 


The total lateral force and center of pressure for the beam 
are calculated by using rigid body translation and rotation modes 
in the following matrix product. 

j^U} 

where 

{1} is a column of 
(PP) is a column of 
{Z} is a column of 

= total force on the beam, and 
P^ = total moment of forces on the beam about x = 0. 


i - ['t 

ones , 

the panel point x stations, 

the concentrated panel point forces. 


From this data the center of pressure of the forces is calculated 
from 


X 


cp 



(6a) 


(6b) 



ALODl— 6/7 


EXAMPLE 

Consider a beam with distributed and concentrated lateral 
forces as shown in the sketch. 



The panel points are denoted by the points numbered 1 thru 5. 
The panel point x-stations and the force data are defir4ed in: 


{PP} 



[DIST] = [l5. 

30. 

5. 

20. 

45. 

55. 

10. 

10. 


and [CONC] = (70, 


30 .]. 


Using the technique described previously and Equations (5) and 
(6), the reader can verify the following results. 
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For row 1 of iDIST]; 

Bay 1; Z] = 27.80 lb 

Z 2 = ''^.20 lb 

Bay 2: z, = 76.07 lb 

Z 3 * 11. A3 lb. 

For row 2 of [DIST]: 

Bay 3: Z 3 = 75.00 lb 

- 25.00 lb. 

For row 1 of (CONC]: 

Bay A: zi, = 20.00 lb 

Z 5 = 10.00 lb. 

The final forces at the selected panel points {PP} to replace 
distributed forces defined by iDlST] and concentrated forces de- 
fined by I CONC] are given by the sum of the above results as 

" 27.80 
1A8.27 
86 . A3 
45. CO 
lU.OO 


Ul = 
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Subroutine AL0D2 Cakes (on option) distributed and concen- 
trated axial forces on a beam and replaces them with representa- 
tive concentrated forces at selected points on the beam. 

The x-stations of the selected points (panel points) are 
given in {PPl. These x-stations must be in increasing order. 

The distributed axial force, f(x), is assumed to be piece- 
wise linear and is represented by straight line segments as shown 
in Figure 1. 

f(x) 



The x-stations of the end points for the line segments giv- 
ing ? ’ke distributed axial force are independent of the panel 
po‘ t x-stations. However, the distributed axial force must be 
within the panel point limits. The line segments representing 
the distributed axial force may or may not be joined and may 
overlap- The distributed axial force is defined in [DIST]. Each 
row of [DIST] represents one nonvertical line segment. The form 
Oi each row of [DIST] is [x] x^ fi ^ 2 ! where xj, fj give the first 
end point and X 2 , ^2 second end point of a line segment. 

concentrated axial forces are defined in [CONC] . Each 

row of [CONC] contains one concentrated axial force, F , and its 

c 

X location in the form jx^ ^ concentrated axial force may 

be outside the panel point limits. 

The calculated representative concentrated forces at the 
selected panel points are placed in {Z}. The total axial force 
is alsj calculated and printed. 
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DESC R IPTION OF TECHNIQUE 

The replacement of distributed and concentrated axial forces 
by representative concentrated forces at selected panel points 
is obtained using a virtual work approach as follows* 


The virtual work done by the axial forces on a beam is de- 
fined by 


AW 



f(x) Ag(x) dx + ^ ^ 
c 


where 

f(x) is the distributed axial force, 

F is a concentrated axial force, 
c 

6g(x) is the axial virtual displacement, 

X is the longitudinal axis of the beam, 
x^ is the starting x-station of the beam, and 
Xg is the ending x-station of the beam. 


( 1 ) 


The finite summation is over the number of concentrated axial 

forces, F . 

c 

Most techniques for describing the axial displacement (vir- 
tual or real) assume a function, g(x), over the entire length of 
the beam. However, considerable skill is required in choosing 
this function. The technique used here is to represent the 
axial displacement between consecutive panel points by some simple 
function of the panel point displacement. A uniform displacement 
function will be assumed here. The displacement between consecu- 
tive panel points k and k+1 is assumed equal to the displacement 
ot panel point k+1. That is. 


g(x) = 8(Vl) ^ 


Similarly 


6g(x) = 6g 6gj^^^. 


( 2 ) 
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The region between panel points k and k+1 is referred to as bay 
k« 


The virtual work of distributed and concentrated forces will 
be considered separately. The distributed axial force, f(x), is 
considered first. The geometry for a line segment of distributed 
force is shown in Figure 2. 


f(x) 



Xi 


Figure 2 Line Segment Geometry 

The equation for a straight line segment as shown in Figure 3 is 

f(x) = fj + (x “ xi) (f^ - fi)/(x 2 - xi). (3) 

Substituting Equations (2) and (3) into (1) gives the virtual 
work of the distributed axial force represented by one line seg- 
ment 1 in bay k as 




^)/{\ - *p)) 


(4) 


The subscripts p and q have been introduced to handle the possi- 
bility of a line segment extending past the bay limits. Thus, 

x is the greater of xi or x, and x is the lesser of or x, . , 
p ® ^ k q 2 k+1 

Similarly, f^ is either f| or f^^ and f^ is either f 2 or 

The integration is continued for the line segment in adjacent 
bays, if necessary, until the entire line segment has been used. 
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Performing the integration of Equation (4) yields 

*"i.k * Vi ‘8k+i 

where 


z, , = * 5 /f + f \ /x - X 

k+1 [p q/ \ q P 


)• 


(5) 


(5a) 


By definition, the coefficient of the virtual displacement is the 
force. Thus, is the representative concentrated force at 

panel point k+1 to replace the distributed axial force represented 
by one line segment i in bay k. 


The concentrated axial forces are considered last. For a 
force located in bay k, the displacement at the force location 
is given by Equation (2) as 



6g 


k+1* 


From Equation (1), the virtual work of one concentrated force, 
F^, in bay k is given as 


c , k k+? 




k+1 


where 


k+1 



'. 6 ) 

(6a) 


The virtual displacement coefficient, is the representative 

concentrated force at panel point k+1 to replace the concentrated 
axial force, F^, in bay k. 

The representative concentrated panel point forces for the 
entire beam are finally obtained by evaluating Equation (5a) for 
each line segment of distributed force and Equation (6a) for each 
concentrated force. Like terms, are then summed. 

It is interesting to note that the virtual work approach used 
here, with displacement in bay k assumed uniform and equal to the 
displacement at panel point k+1, gives identical results to the 
more common procedure of placing all forces in bay k at panel 
point k+1. 

The total axial force is calculated by summing the elements 
of (Z). 
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EXAMPLE 


Consider a beam with distributed and concentrated axial forces 
as shown in the sketch. 



I 1 1 > 1 I I 1 i I 1 11 -i- 1- 

0 10 20 30 40 50 60 70 


X. 

80 


X (in.) 


The panel points are denoted by the points numbered 1 thru 5. 
The panel point x-statlons and the force data are defined In: 


{PP} 


10 . 

25. 

45. 

65. 

80. 

_ _ 


[DIST] = 


15. 30. 5. 20. 
45. 55. 5. 5. 


and [CONC] = [70. 


30.]. 


Using the technique described previously and Equations (5a) 
and ( 6 a), the reader can verify the following results. 


For row 1 of [DIST]: 

Bay 1: Z 2 = 100.0 lb 

Bay 2: Z 3 = 87.5 lb. 

For row 2 of [DIST]: 

Bay 3: zt, = 50 lb. 

For row 1 of [CONC]: 

Bay 4: 25 = 30.0 lb. 
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The final forces at the selected points (PP) to replace distributed 
forces defined by [DIST] and concentrated forces defined by [CONC] 
are given by the sum of the above results as 


100.00 

87.5 

50.0 

30.0 
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ASSEM 


Subroutine ASSEM places (assembles) a matrix [A] into a sec- 
ond matrix [Z] starting at a designated row, column location 
(IRZ, JCZ, respectively) in [Z]. The elements of [A] will replace 
corresponding elements in the original [ZJ. Before the first use 
of this subroutine in forming IZ), it is important that [Z] is 
correctly defined. For example, if [Z] is to be originally all 
zeros* Subroutine ZERO could be used. This subroutine may be 
called repeatedly to form [Z] from the assembly of several [A] 
matrices. The [A] matrix must be within the row, column limits 
of [ZJ. In subscript notation, 

/k = 1. NRA\ 

ij k«. \«- = 1, NCA/ 


where 


1 - k + IRZ - 1 
j = + JCZ - 1 

NRA Is the number of rows of [A] ; 
NCA is the number of columns of [A] . 
EXAMPLE 

Consider a matrix defined as 






1. 

0. 

0. 

0. 



[Z] 

3x4 

s: 

0. 

2. 

0. 

0. 






0. 

0. 

0. 

0. 



("3. 

4. 


5.1 





Matrix 

■ [e. 

7. 


8.J 

is 

to be 

assembled into [Z] start 

Ing at 

the 2,1 location 

(IRZ - 

2. 

JCZ » 

1) of [Z]. 


The result of this operation will be 



1 . 

0. 

0. 

0. 

II 

3. 

4. 

5. 

0. 


6 . 

7. 

8 . 

0 . 



ATXBAl 

Subroutine ATXBAl calculates a special matrix product. In matrix 
notation, 

T 

^^^NCBxNCB ^NCBxNRB ^®Wb>'NCB 

where [Z] is calculated as symmetric and placed in the same core 
locations as [A] to allow larger matrix sizes. NRB is the number 

T 

of rows of [B], columns of [A] . NCB is the number of c^jLumnb 

T 

of [B], rows of [A] , and size of [Z]. 

In scalar summation notation, 

NRB NRB 

^Ij " ^ik ‘^kj “ S ^ki 

k=l k=l 

DESCRIPTION OF TECHNIQUE 

To reduce computer time and accomplish the matrix product using 
only two matrix core spaces, an intermediate work space vector 
is used. The size of this work vector determines the maximum 
number of columns of [B], i.e., NCB. The matrix multiplication 
is accomplished as follows. A single row of [A]^ (columns of [A]) 
is multiplied times the columns of [B]. These results are stored 
in the work vector until all the columns of [B] have been used. 

Because [Z] will be symmetric, columns of [B] less than the row of [A] 
need not be used. The elements in the work vector then replace the 
row of [A]^ (column of [A)) used in the multiplication. This 
procedure is repeated for all rows of [A]^ (columns of [A]) to 
calculate the lower half of the answer [Z]. The lower half is then 
reflected to the upper half to give the final result. 


/i = 1, NCB\ 
\j - I, NCB/ 



ATXBB 


Subroutine ATXBB calculates a special matrix product. In matrix 
notation, 

T 

^^^NRATxNCB ^ ^NRATxNRB "^^NR3>:NCB 


where [Z] is placed in the same core locations as [B] to allow 

T 

larger matri:: sizes; NRAT is '*e number of rows of [A] and [Z]; 

T 

NRB is tho. number of rows of [B] and columns of [A] ; NCB is the 
number of columns of [B] and [Z]. 

In scalar .ummation notation 

1, NRAT 
1, NCB 

DESCRIPTION OF TECHNIQUE 


'ij 


NRP 

E 

k=l 




NRB 


k=l 


^ki ^kj 




To reduce computer time and accomplish the matrix product using 
only two matrix core spaces, an intermediate work space vector 
is used. The size of this work vector determines the maximum 

T 

number of rows of [A] > i.e., NRAT. The matrix multiplication is 

.T 

accomplished as follows. The rows of [Aj (columns of [A]) are 
multiplied times a single column of [B]. These results are stored 

T 

in the work vector until all the rows of [A] have been used. The 
elements in the work vector then replace the column of [B] used in 
the multiplication. This procedure is repeated for all columns 
of [B] to calculate the result [Z]. 



Subroutine ATXBBl calculates a special matrix product 
matrix notation. 


In 


ATXBBl 


tz] 


NRBxNCB 


= (lA]’') 


NRBxNRB 


[B] 


NRBxNCB 


where 


[A] = a square, upper triangular matrix; 

[Z] = placed in the same core locations as [B] to allow larger 
m :rix sizes; 

NRB = the number of rows of [B] and [Z] and size of [A]; 

NCB s the number of columns of [B] and [Z]. 

In scalar summation notation 


z 


ij 


NRB 


NRB 

E 

cr 

II 

E 

k«l 


k=l 


a, . b, . 
ki kj 


/i = 1, NRB\ 

\j = 1, ncb/ 


DESCRIPTION OF TECHNIQUE 

To reduce computer time and accomplish the matrix product using 
only two matrix cere spaces, an intermediate work space vector 
is used. The size of this work vector determines the maximum 
NRB. The matrix multiplication is accomplished as follows. The 
T 

rows of [A] (columns of [A]) up to the diagonal are multiplied 
times a single column of [B[. These results are stored in the 

T 

work vector until all the rows of [A] have been used. The ele- 
ments in the work vector then replace the column of (Bj used in 
the multiplication. This procedure is repeated for all columns 
of [B] to calculate the result [Z]. 
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Subroutine ATXBB2 calculates a special matrix product. In 
matrix notation, 

T 

^^^NCBxNCB ^ ^NCB>‘NRB ^^^NRB’^NCB 

where 

[A[ - rectangular matrix; 

[B] > rectangular matrix; 

[Z] = square, symmetric matrix; 

NRB > number of rows of [B] and [A]; 

NCB s number of columns of [B] and [A]. Number of rows and 
columns of [Z]. 

The result, [Z], is formed in the same core location as [B]. 

The algorithm for the elements of the matrix product, [Z], is 
NRB NRB 

‘ij * "ji ■ L \i - E V '■kj (j ; i: SS) 

k=l k=l 

DESCRIPTION OF TECHNIQUE 

This subroutine is Intended to be used as the second multiplication 
of a triple matrix product 

[Z] = [A]*^ [C] [A] 

and matrix [B] represents the product of [C] and [A] 

[BJ = IC] [A]. 

Matrix [C] is assumed to be symmetric and it follows, therefore, 
that [zj is also symmetric. The upper half of [Z] is calculated 
and reflected to form the lower half. The result, [Z], is placed 
in the [B] core locations. 
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SubrouCine AXBAl calculates a special matrix product. In matrix 
notation, 

^^^NRAxNCA “ ^®^NCAxNCA 

where 


[A] ” rectangular matrix; 

(B] •= upper triangular, squ re matrix; 

[Z] rectangular matrix; 

NRA » number of rows of [A] and [Z]; 

NCA number cf columns of [A] and [Z] ; 

number of rows and columns of [B]. 

The result, [Z], is formed in the same core locations as [A]. 

The algorithm for the elements of the matrix product, [Z], is 


i 



k-1 


/i - i, NRA\ 
\j = 1, NCA/ 





AXBA2 


Subroutine AXBA2 calculates a special matrix product, in matrix 
notation, 

where 

[A] » square matrix; 

[B] « upper triangular, square matrix; 

[Z] a symetrlc, square matrix; 

N * number of rows and columns In [A], [B], and [Z]. 

The result, [Z], Is formed In the same core locations as [A]. 

The algorithm for the elements of the matrix product, [Z], Is 

(j-i:") 

k=l 



Subroutine AXBA3 calculates a special matrix product. In matrix 
notation. 


AXBA3 


^^^NRBxNRB ^®^NRB^NCB 


where 

[A] = upper triangular, square matrix; 

IB] = rectangular matrix; 

[Z] = rectangular matrix; 

NRB - number of rows of [B] and [Z], 

number of rows and columns of [A] ; 

NCB « number of columns of [B] and [Z]. 

The result, IZ], is formed in the same core locations as 

The algorithm for the elements of the matrix product, [Z], is 


*ij 


NRB 



k=i 


/i = 1, NRB\ 
\j = 1, NCB/ 


*ik '’kj 
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Subroutine BABT calculates a special form of the triple matrix 
product. In matrix notation. 


^^^NRBxNRB “ ^°^NRBxNCU ^^^NCBxNCB ^^^NCBxNRU 


where 


z 


ij 


MCB NCB 

~ ^ ^it 

t=l k=l 


/i = 1, NRB\ 

\j = 1, nrb/ 


CAj » square, non- symmetric matrix. 

Cb 7 * rectangular matrix. 

(Zj » square, non- symmetric matrix. 

NRB = number of ro’vs of C®] and size of [zJ. 

NCB » number of columns of [B] and size of [aJ. 

Theorem ; If [A] is symmetric jthat is, [A] = [A]^) then [Z] is 
symmetric. 

X 

Proof ; [Z]^ = ([B] [A] t®]'*’) 

= lA]’^ 

= IBI lA] 

= tz]. 

Therefore, if [Aj is non-symmetric, [z] is non- symme trie . 
DESCRIPTION OF TECHNIQUE 

To reduce computer time, an intermediate work space vector is 
used in the subroutine. The size of this work vector determines 
the limitation on the number of columns of [B] and thus the size 
of (A] (i.e., NCB). This special triple matrix product is accom- 
plished as follows. The rows of (Aj are mulLlplled times a 

column of [Bj^. These results arc stored in the work vector. 

Next, the rows of [b] multiply the work vector to obtain a column 
of Che answer fZj. This procedure is reported for all the columns 
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T 

of [B] to give the answer [z3 . 

EXAMPLE 

Consider input of 

7. 10. 

- 8 . 11 . 

9. 12. 

The reader can easily verify the output to be 


CA] 


2x2 


13. 

14. 

15. 

16. 


and (eQ 


3x2 


[Zl 


3*3 


4267. 

910. 

5265 

1067. 

216. 

1317 

5259. 

1122. 

6489 
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Subroutine BABTA calculotes a special form of the triple matrix 
product. This subroutine is a modification of Subroutine BABT to al> 
low larger matrix sizes by placing the answer [Z] in the same core 
locations as [A] . In matrix notation, 


^^^NRBxNRB ” ^^^NRBxNCB ^^^NCBxNCB ^‘^^NCBxNRB 


where 


'ij 


MCB NCB 

^ii ^£k ^jk 

t=l k=l 


/i = 1, nrb\ 
u - 1, nrb/ 


[a] = square, non- symmetric matrix. 

[B] > rectangular matrix. 

[Z] B square, non- symmetric matrix. 

NRB number of rows of [b] and size of [z]. 

NCB » number of columns of [B] and size of [A). 

Theorem ; If [A] is symmetric [that is, [A] = [A]^) , then [Z] is 

synatetric. 

T 

Proof ; IZ]^ = ([B] [A] [B]'^') 

» ([B]^V [A]^ [B]"^ 

“ [B] [A] [B]^ 

= [ 2 ], 

Therefore, if [A] is non- syirane trie , [,Z] is non-symmetric . 
DESCRIPTION OF TECHNIQUE 

To reduce computer time and to accomplish tlie special triple matrix 
product using only two matrix core spaces, an intermediate work space 
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vector is used in the subroutine. The size of this work vector deter- 
mines the limitation on the number of columns of [B] (i.e., NCB) . Tlte 
special triple matrix product is accomplished as follows. A row of [Aj 
is placed in the work vector whicl) is multiplied times the columns of 
[B]T. These results then replace the row of [a] used in the multipli- 
cation. This procedure is repeated for all rows of [A] to obtain [AB^] . 
Next, a column of [AB^] is placed in the work vector which is multiplied 
by the rows of fBj to tilve a column of the answer [Z) . This last pro- 
cedure is repeated for all columns of [AB^) to give the answer fz] . 

EXAMPLE 

Consider input of 

14. 

X.U » 

The reader can easily verify the output to be 

4267. 910- 5265. 

[Zl, , = 1067. 216. 1317. 

3x3 

5259. 1122- 6489. 






BTAB --1/2 


Subroutine BTAB calculates a special triple matrix product. In 
matrix notation, 


^^^NCBxNCB 


[B] 


T 

NCBxNRfi 


^^^NKBxNRB ^®^NRBxNCB 


where 


'ij 


NRB NRB 

^Hi *^)ik ^kj 

«.»1 k=l 


/ i - 1 , NCB \ 

U = 1, ncb/ 


[a] ■ square, non- symmetric matrix. 

[B] « rectangular matrix 

[Z] ® square, non- symne trie matrix 

NRB = number of rows of Qs] and size of [A] . 

NCB = number of columns of [B] and size of [z] . 

Theorem ; If [A] is symmetric (that is, [A] = [A]^ ), then [Z] is 
symmetric . 

[Z]^ = ([B]^ [A] [B])^ 

= [B]" lA]" ([B]'^)" 

= [B]'^ [A] [B] 

= [Z]. 

Therefore, if [A] is non- symmetric, [z] is non-symuetric . 
DESCRIPTION OF TECHNIQUE 


To reduce computer time, an intermediate work space vector is 
used in the subroutine. The size of this work vector determines the 
limitation on the number of rows of [B] and thus the size of [A] (i.e., 
NRB). The triple matrix product is accomplished as follows. The rows 
of [AJ are multiplied times a column of [b] . These results are stored 



BTAB -2/2 


in the work vector. Next, the rows of (columns of (bJ) multiply 

the work vector to obtain a column of the answer [z] . This procedure 

is repeated for all the columns of M to give the answer [z] .. 


EXAMPLE 

Consider Input of 



1 — ~ 

>— < 

1 


[7. - 

-8. 9.1 

^^^2*2 

15. 16. 

- - 

and [B]2^3 = 

10. 11. 12. 

The reader can easily verify the output to be 




*4267. 910. 

5265.' 




1067. 216. 

1317. 




5259. 1122. 

6489. 

• 
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Subroutine BTABA calculates a special triple matrix product. This 
subroutine is a modification of Subroutine BTAB to allow larger nuitrix 
sizes by placing the answer [2j in the same core locations as [a] . In 
matrix notation. 


^^^NCBxNCB “ ^^^NCBxNRB ^''^NKUxNRB ^“^NRBxNCU 


where 


z 


ij 


NRB NRB 


EE 


b . 


Ilk 



/i = 1, NCb\ 
\j = 1. NCB/ 


M " Square, non- symmetric matrix. 

[b] = rectangular matrix. 

Cz] “ square, non- symmetric matrix. 

NRB « number of rows of [B] and size of (VQ • 

NCB *> number of columns of [Bj and size of [z] . 

Theorem ; If [A] is symmetric (that is, [A] = [A]^), then [Z] is 
symmetric . 

IZ]'^ = (iBj'^' [A] [B])^ 

T 

= [A]'^ ([B)M* 

= iB]^ lA] [B] 

» [ZJ. 

Therefore, if [aJ is non- symae trie, [z3 is non- symme trie . 
DESCRIPTION OF TECHtUQUE 

To reduce computer time and to accomplish the triple matrix product 
using only two matrix core spaces, an intermediate work space vector is 
used in the subroutine. The size of this work vector determines the 
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limitation on the number of rows and columns of [uj (i.e., NltB and NCB) 
The triple matrix product is accomplished as follows. A row of [A, is 
placed in the work vector which is multiplied the columns of [B^ 

These results then replace the row of [Aj used in the multiplication. 
This procedure is repeated for all rows of [a] to obtain [AB^ . Next 
the rows of (columns of (bJ) are multiplied times a single column 

of [AB] just formed. These results are stored in the work vector until 
all the rows of [B]*^ have been used. The elements in the work vector 
then replace the column o^ [AB] used in the multiplication. This pro* 
cedure is repeated for all columns of [AB] to give the answer [Z, .. 


LX/V'IPLL 


Consider input of 




- 


r‘ 



13. 

14. 


7. 

-8. 9. 

■ 

15. 

16. 

■J 

and [ 012 x 3 “ 

_i0. 

11. 12._ 


The reader can easily verify the output to be 


4267. 

910. 

5265. 

1067. 

216. 

1317. 

5259. 

1122. 

6489. 
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Subroutine BTABA2 calculates a soecial tripia matrix product. In 
matrix notation. 


IZ] 


NxN 




[B] 


NxN 


where 

(A] = square, symmetric matrix; 

[B] = square, upper triangular matrix; 

[Z] = square, symmetric matrxx. This answer is placed in the 
same core locations as [A] to allow larger sizes; 

N * size of each matrix, 

T 

Theorem: T.f [A] is symmetric (that is, [A] = [A] ), then [Z] is 

synonetric. 

Proof: 

tzf = [A] [B])^ 

- [B]' [A]^ ([B]") 

* [B]'^ [A] [B] 

= [Z]. 

DESCRIPTION OF T ECHNIQUE 


To reduce computer time and accomplish the triple matrix product 
using only two matrix core spaces, an intermediate work space 
vector is used. The size of this work vector determines the 
maximum matrix size, i.e., N. The triple matrix product is 
accomplished as follows. A single row of [A] is placed in the 
work vector and is multiplied times the columns of fb]. These 
results then replace the row of [A] used in the multiplication. 
This procedure is repeated for all the rows of [A]. Use is made 

T 

of the upper triangular form of [B]. Next, the rows of [B] 
(columns of [B]) are multiplied times a single column of [AB] 
just formed. These results are stored in the work vector until 

T 

all of the rows of [B] have been used. Again use is made of the 
upper triangular form of [B]. The elements in the work vector 
then replace the cclumn of [AB] used in the multiplication, down 
to the diagonal and the corresponding row out to the diagonal. 
This procedure is repeated for all columns and corresponding rows 
of [AB] to give the answer [Z]. 
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Example: 

Consider input of 



The reader can easily verify the result to be 




COLMLT 


Subroutine COLMLT evaluates a special matrix operation by mul- 
tiplying each column of a matrix [fi] by a scalar. That is, 

^^^NRxKC “ p^^^^^NRxi 

where 


z . . 
^2 


a. b . , 




denotes column j of [B] . 


Each scalar 


is an element of 


the input vector {AVECi. NR is the number of rows in IB] and 
[Z], and NC is the nunber of columns in [B] and [Z] and the size 
of (AVEC). The number of elements of {AVEC} must be equal to the 
number of solumns of [BJ* 


EXAMPLE 


Consider input of 

(AVEC) = [2. -3. 4J and 



The output will then be 






COMENT 


Subroutine COMENT reads input comment cards and reproduces 
each card in the printed output of the computer run* Each com- 
ment card may have any keypunch symbol in card columns 1 thru 
78. A use of COMENT is to print an explanation of coordinates 
used in a computer run* Thus, this information is always re- 
tained with a run to correlate matrix location numbers with 
physical coordinates* 
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Subroutine DCOMl decomposes [A] to form an upper triangular 

T 

matrix [Z] such that [A] = [Z] [Z]. [A] must be real, square, 

symmetric and positive definite* The Choleski square root method 
is used* In matrix notation. 




where 


^^^NxN ~ 

’^11 

^12 

••• 


0 

^22 

• ' * ^2H 


0 

0 

• • • 


N is the size of the matrices (square) . 


DESCRIPTION OF TECHNIQUE 

To determine the elements of [Z], consider the multiplica- 
tion of two matrices* 

[A] = [B] [Z], 

A general term of [A] is 

•ij ■ '>11 hi * ‘>12 ^2j ■" ••• ■"'“in "nj- 


But 

[B] = [Z]^, 


thus 


b 


ij 


z 


ji* 


Therefore 


^2i ^2j "^ * * * ■'■ ^Ni ^Nj • 



DCOMl—2/3 


Because 


ve have 


'ij 


0 for i ^ j , 


and 


a. . = z- . 2, , + . + 

ij Ij 2i 2j 


+ z < 


IX ij 


a. . 
11 




+ 


+ z 


2 

ii* 


(i - j) 


From these last two equations the formulas for determining z.. 
are obtained. That is, 



MISCELLANEOUS 

The diagonal elements of [Z] are the square root of the deter- 
minant ratios of [A] . The determinant of [A] is the product of 
these determinant ratios. That is, 

N 

iAi - n 

i=i 
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EXAMPLE 

The input matrix to be decomposed is 


[A] 


3x3 



Using the technique described previously, the reader can 
verify the output to be 



REFERENCE 

Faddeeva, V. N. : Computational '-lethods of Linear Algebra. 

Dover Publications Inc., New York, 1959. 



Subroutine DI^G places the elements from a vector (row or col- 
umn matrix) |AVEC| into the corresponding diagonal locations of a 
square matrix [Z] and sets the off-diagonal elements of the square 
matrix to zero- In subscript notation. 


11 ■ “i 

(i 

= 1, N) 

IJ ■ 

(1 

* j) 


In matrix notation. 


[Z] 


NxN 





where N is the size of [Zj (square), and the length of ^AVEC 
EXAMPLE 


Consider input of 


|avec| = 


1 . 
2 . 
3. 

and N = 3 . 


[Z], 


3x3 


will 

give 


1. 

0. 

0, 

0. 

2. 

0. 

0. 

0. 

3. 


U1A( 
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Subroutine DIFFl performs numerical differentiation assuming 
a linear function between known points. The x and y coordinates 
of the known points are given by the elements of {XA} and the cor- 
responding elements in a column of [YA], respectively. Kach column 
of [YA] gives the y coordinates of a different set of points. 
Derivatives are calculated at selected x coordinates whidi are 
given by the elements of {XZ}. These derivatives are placed in 
[Z]. Each column of [Z] has derivatives of the respective column 
of [YA]. Differentiation of an extrapolated linear function is 
performed when any element of {XZ} exceeds the limits of fXA). 

DERIVATION OF TECHNIQUE 

Giveti the x,y coordinates of points i and i+1 , the derivative 
dy/dx at X, is to be found by assuming a linear function. 


y 



From the above sketch it is obvious that 


dx 



( 1 ) 


The following tabulation gives the correlation between the nomen- 
clature of Equation (1) and that used in the Fortran coding in 
the subroutine. 


Equation (1) 

Fortran Coding* 


XA(I) 

’^i+i 

XA(I+1) 


YA(I) * 


YA(I+1) * 


XZ (K) 

dy/dx 1 

Z(K) * 

■ 



’^SuhscrlpL J, dcntUlng dlffcront sets of points, 
lias been omitted iur clarity. 
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EXAMPLE 


Consider the folloii^ing two sets of points denoted by (1) and 


( 2 ). 



The coordinates of the points are given by 



-2. 


-A. 

2. 

{XA} = 

1. 

and [YA] = 

2. 

-4. 


3. 


1. 

-2. 


{XA} gives the x coordinates of the points in both sets (1) and 
(2). Column 1 of [YA] gives the y coordinates of the points in 
set (1) and column 2 of [YA] gives the y coordinates of the ooints 
in set (2). Derivatives are wanted at x = -1. and x = 7., that is. 


at {XZ} = 


- 1 . 

7. 


At X = ~1*» derivatives of dy/dx = 2. and 


dy/dx = -2. are calculated from columns 1 and 2 of [YA], respec- 
tively, using rows 1 and 2 of {XA} and [YA]. At x = 7., deriva- 
tives of dy/dx = -0.5 and dy/dx = 1. are calculated from columns 
1 and 2 of [YA], respectively, using rows 2 and 3 of {XA} and 
[YA]. The final result is [Z] = 


2 . 

-0.5 


-2. 

1 . 


To relate this example problem to a practical problem, con- 
sider (XA) to be the collocation points (panel points) of a vehicle 
and [YA] to be the modal displacements for two modes. Gyros are 
to be placed at stations x = -1. and 7. (i.e., {XZ}). The modal 
slopes ([Z]) are to be found at these stations. 
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Subroutine DIFF2 perforius numerical differentiation assuming 
a diparabolic function between known points. A parabolic function 
is used where only three points are av^-ilable. A diparabolic 
function is obtained from the weighted average of two adjacent 
parabolas and will be explained later. The x and y coordinates 
of the known points are given by the elements of {XA} and t!ie cor-* 
responding elements in a column of [YA], respectively. Fach column 
of [YA] gives the y coordinates of a different set of points. 
Derivatives ar^ calculated at selected x coordinates which are 
given by the elements of IXZ}. These derivatives are placed in 
[Z]. Each column of [Z] has derivatives of the respective column 
of [YA]. Differentiation of an extrapolated parabolic function 
is performed when any element of {XZ} exceeds the limits of (XA). 

DERIVATION OF TECHNIQUE 

The diparabolic differentiation procedure is obtained as fol- 
lows. Because this procedure is dependent upon using parabolas, 
the parabola will be considered first. Given the x,y coordinates 
of points 1, 2, and 3, a parabola is to be fitted to tliese points. 



The equation for a parabola wi;h axis parallel to the y axis is 

y (x) = Ax^ + Bx + C 


or 


where 


or 


y(H) = [h2 H 1] 


a 

b 

c 



(la) 


(lb) 
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is used for ease in later algebraic calculations. The coefficients 
a, b, and c can be determined because x (or H) and y coordinates 
at poiats 1, 2, and 3 are known. That is, 



H2(H2-H3) - H|(Hi-H3) + H§(H3- H^) 


therefore 

y(H) = [r 2 H 1] [ij;] 


For a given set of points as shown below, a parabolic function 
is used to the left of point 1 and between points 1 and 2. Also 
a parabolic function is used to the right of the last point n 
and between points n-1 and n. Diparabolic functions are used 
between all other points. 





*k 


The 'derivative dy/dx at is calculated using the chain rule for 
differentiation. That is, 

dy _ dy dH 
dx dH dx 

From Equation (3), 


dH 
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and from Equation' (la). 


- dH 


dx “ Xj 


therefore 


=: 1 
dx X 2 “ Xj 


[2H 1 0] M 


Also, from Equation (la) > 

*l)/(*2 - *1) 

“ (*1 - *l)/(*2 - ’‘1) “ ° 
H 2 - iJXj - *l)/(x2 - *l) = 1 


H 


3 ‘ (’'3 ■ *l)/f 


X2 - *1) "13. 


(4) 


(4a) 


(4b) 


Using these expressions with Equations (Xi and (4), the final 
result is obtained as 


dx 


■= ^ 4 

1/D 

1/(1-D) 

-1/D(1-D) 

X2 - XI 1 k 1 

\ 

-d+D)/D 

-D/(l-D) 

1/D(1-D). 


''1 

V2 

Y3 


(4c) 


The following table gives the correlation between the nomen- 
clature of Equations b> c) and that used in the Fortran coding 

in the subroutine. 


Equations (4a,b,c) 

Fortran Coding’^ 

X 

m 


XA(ra) 



YA(m) * 

\ 


H 



XZ(K) 

dx 


Z(K) * 


*Subscript j, denoting different sets of 
points, has been omitted for clarity. 
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The derivative dy/dx at is found as follo\s. A diparabolic 

function in bay i-1 of the above sketch is obtained as the weighted 
average of parabolas A and B. That is 


y 

(H) = (1-H) + Hy^ 

(5) 

where 

X - X. , 

H 

X, - X. - 
1 1-1 

(6a) 

as in either Equations 

(la) or (lb). 


For parabola A: * 

(“ 1-2 - *i-i)/(‘i - * 1 - 1 ) - ^ 

(6b) 

-1 “ 

(‘ 1-2 ■ ’‘i-i)/(‘‘i ■ ‘i-i) ■ “ 



“l • fl - -l-l)/!*! - ’'l-l) ■ ^ 

From Equations (2) and (3), 

"^1-2" 

y^(H) = (h2 h 1] 

_ — 
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where 



-1/C(1-C) 

1/C 

1/(1-C) 

m 

1/C(1-C) 

-d+O/c 

-C/(l-C) 


0 

1 

0 


For parabola B: ■ ^*^.1 “ *i-l)/^^*i ” ^i-l) " ^ 

"i * ('i - Vi)/(*t ■ -i-i) * * 

"i+i * (*1+1 " *i-i)/(''i “ *i-i) ‘ 

From Equations (2) and (3) » 

'yi-i' 

yg(H) = [H^ H 1] 

7i+l_ 

where 



l/D 

1/d-D) 

-1/D(1-D) 

x; 

-d+D)/D 

-D/d-D) 

1/D(1-D) 


1 

0 

0 


Substituting these expressions for and into Equation (5) 
gives 


= [H^ H 1] Id*] 
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where 




1 

C - D 

D - C 

-1 

C(l-C) 

CD 

(l-O(l-D) 

D(l-D) 

-2 

2D - C 

1-2D + C 

1 

C(l-C) 

CD 

Cl-O(l-D) 

D(l-D) 

1 

-d+c) 

-C 


C(l-C) 

C 

1-C 

U 

0 

1 

0 

0 


The derivative dy/dx at is calculated using the chain rule for 
differentiation* That is» 


dy _ ^ ^ 
dx dH dx 


From Equation (7), 


^i-2 


dH 


= [3H^ 2H 1 0] [ij)] 


^i-l 




^i+1 


and from Equation (6a) 

dH 1 


dx X. - X. - 
1 1-1 


/•At X 




y 


i-2 


iJL 

dx 





(a; 


where Is given above. 
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The following table gives the correlation between the nomen^ 
clature of Equations (6a, b, c and 8) and that used in tlie Fortran 
coding in the subroutine. 


m = 2, 

m * 2, U 0, -I 


^Subscript denoting different sets of points, 
has been omitted for clarity. 



\ 


The derivative dy/dx at is calculated using the chain rule for 
differentiation. That is, 

dy ^ dy dH 
dx dH dx 



From Equation (3) 
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and from Equation (lb) , 


in = 1 

dx X - X , 
n n-1 


therefore 


= 1 

dx ~ X - X , 
n n-1 


[2H 1 0] [.;■] 


^n-2 


^n-1 


(9) 


Also, from Equation (lb) , 

"k - - Viyl-n - Vl) 

“n-2 ' (V2 - Vl)/(’n " Vl) ' 
“n-1 ■ (Vl ■ Vl)/(‘n ■ Vl) ■ “ 
*n • (*« - VlJ/f-n - Vl) ■ 1 • 


Using these expressions with Equations (2) and (9), the final 
result is obtained as 


(9a) 

(9b) 


dx 


1 = — ^ — f^Hi, d 

-1/C(1-C) 

1/C 

l/(l-c) 

|«^ \ ■ *n-l 1 “ J 

1/C(1-C) 

-d+o/c 

-C/(l-C) 


^n-2 


^n-1 


(9c) 


The following table gives the correlation between the nomen- 
clature of Equations (9a, b, c) and that used in the Fortran coding 
in the subroutine. 


Equations (9a,b ,c) 

Fortran Coding* 

X 

XA(NXA-m) 

n-m 

y « 

YA(NXA-m) * 

n-m 

H, 

H 

k 



XZ(K) 

^1 

dx 1 

Z(K) * 

\\ 



*Subscript i, denoting different sets of 
points, has been omitted for clarity* 
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EXAMPLE 


Consider the following two sets of points denoted by (1) ana 



10 r 


1. -1.*” 

20. 


2. -2. 

50. 

and [YA] = 

3. -6. 

70. 


9. -2. 

110^ 


_5. 2^ 


{XA} gives the x coordinates of the points in both sets (1) and 
(2). Column 1 of [YA] gives the y coordinates of the points in 
set (1) and column 2 of [YA] gives the y coordinates of the points 
in set (2). Derivatives are wanted at x = x = 100., and 
X * 65. , that is, at 


{XZ} 


5. 

100 . 

65. 


At X - 5., derivatives of dy/dx = 0.1333 and dy/dx = -0.08333 are 
calculated from columns 1 and 2 of [YA], respectively, using rows 
1, 2, and 3 of {XA} and [YA]. At x = 100., derivatives of dy/dx 
= -0.2333 and dy/dx = 0.0666 are calculated from columns 1 and 2 
of [YA], respectively, using rows 3, 4, and 5 of (XA) and [YA]. 

At X = 65., derivatives of dy/dx = 0.3083 and dy/dx = 0.2354 are 
calculated from columns 1 and 2 of [YA] , respectively, using 
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rows 2, 3, 4, 


and 5 of 


tzj = 


{XA} and [YA]. The result is 

0.1333 -0.08333' 

-0.2333 0.0666 

0.3083 0.2354 


To relate this example problem to a practical problem, con- 
sider {XA} to be the collocation points (panel points) of a 
vehicle and [YA] to be the modal displacements for two modes. 
Gyros are to be placed at stations x = 5., 100., and 65. (i.e., 
{XZ}). The modal slopes ([Z]) are to be found at these stations. 


REFERENCES 

Griffin, J.A. i "A Diparabolic Method of Four-Point Interpola- 
tion." Journal of the Aeronautical Sciences ^ Vol. 28, No. 2, 
Readers' Forum, February 1961. 
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Subroutine DISA removes (disassembles) a matrix [Z] from 
matrix [A] starting at a designed row, column location (IRA, JCA, 
respectively) in LA]. The IZ] matrix must be within the row, 
column limits of (Aj. In subscript notation, 

/i = 1, NRZ\ 

^ij \j = 1, NCZ/ 


where 


k » i + IRA - 1 
£ = j + JCA - 1 

NRZ is the number of rov;s of [Z], and 
NCZ is the number of columns of [Z]. 

EXAMPLE 

Consider a matrix defined as 

1 . 0 . 0 . 0 . 

3. 4. j. 0. 

6. 7. 8. 0. 

Matrix to obtained from [A] starting at the 2,1 loca- 

tion (IRA = 2, JCA r 1) of [A]. The result of this operation will 
be 

3. 4. 5. 

6. 7. 8. 

The matrix [A] remains as originally defined. 
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Subroutine EIGNl calculates the eigenvalues (proper values, 
characteristic roots, latent roots) and eigenvectors (proper 
vectors, characteristic vectors, latent vectors) of a real sym- 
metric matrix using a method of C. G. J. Jacobi. The eigenvalue 
problem can be expressed as 

lA] [^] =1^] t A 3 (1) 


where 


[A] is a real symmetric matrix of order N, 

[4»] is a matrix whose columns are the eigenvectors of 
Equation (1) , 

[ A 3 is a diagonal matrix whose elements are the corre- 
sponding eigenvalues. 

This method diagonalizes matrix [A] by successive plane rota- 
tions and leads simultaneously to all eigenvalues and eigenvectors. 
Theoretically, an infinite number of plane rotations are necessary 
to produce the diagonal form with all off-diagonal elements equal 
to zero. Practically the number of plane rotations is limited to 
a finite number by stopping the process wheii the off-diagonal ele- 
ments are less than a prescribed value. This prescribed value 
(denoted as FOD in the subroutine) may be selected and input by 
the analyst or at the analyst's option calculated in the sub- 
routine as a constant times the trace (sum of the diagonal ele- 
ments) of [A]. 

A threshold version is used. That is, eacih off-diagonal ele- 
ment of [A] is compared in regular sequence with a threshold 
value and a plane rotation performed only if its magnitude ex- 
ceeds the threshold value. This is done to give faster conver- 
gence to the diagonal form. Since [A] is symmetric only the upper 
half need be examined. The threshold value is lowered whenever 
there are no remaining off-diagonal elements with magnitudes 
larger than the threshold. The lowering of the threshold is con- 
tinued until the threshold is less than a prescribed value (FOD) • 

The initial threshold value in the subroutine is obtained by 
dividing the maximum off-diagonal element o.: [A] by 10- Sub- 
sequent reductions in the threshold value aie accomplished by 
dividing the threshold value 10 also* 
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The threshold version is the fastest version of Jacobi's 
method (Reference 2) . 

An important feature of the Jacobi method is that equal 
eigenvalues and the corresponding eigenvectors are easily cal- 
culated with no change in the algoritlimv 

DESCRIPTION QF TECHNIQUE 

Assume that an N^^ order matrix (T] exists that is ortho- 
normal that is,([T]^ [T] = [I]j such that 

[T]^ [A] [T] = t D J, (2) 

where [ i is a diagonal matrix and [A] is the original matrix* 
Premultiply Equation (2) by [T] to give 

[A] [T] = [T] r D 4 . (3) 

Comparing Equations (1) and (3) it is seen that if the matrix [T] 
can be found which transforms a real symmetric matrix [A] into a 

diagonal matrix t ^ then the i^^ diagonal element of D J is 

th th th 

the i eigenvalue of [A], and the i column of [T] is the i 

eigenvector of [A] . 

The orthonorraal matrix [T] is built in stepwise fashion using 
elementary orthonormal transformations 1^]^^ to annihilate, in turn, 

Sr.l^cted off-diagonal elements of [A]* That is, starting with the 
;^iven matrix [Aj , operate as follows, 

[A] = [A]q 

[T]J [A]q [I]^ = (A]^ 

[T]2 [A]^ IT]2 = Ul2 



[A] 


k-1 


[T], 


[A], 
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so that 

[T]j^ . . . [T]2 [T]^' [A] [Tl^ IT]^ [T]j^ = [A]j^ . (4) 

Comparing Equation (4) with Equation (2), the must be chosen 

so that the triple matrix product result [A]^ converges to 

f D ] =* t ^ J ^nd the continued product [T]^ ^^^2 

verges to [T] = [4»]. The most important property of [T]. is that 

K 

each triple matrix product (rotation) with in Equation (4) 

causes a particular off-diagonal element in [A]^ ^ to vanish. As 
mentioned before, an off-diagonal element of lA]^ ^ is picked for 

elimination if it exceeds a threshold value and is referred to as 
the pivot element* If the p, q element of lA]j^ ^ meets the cri- 
teria, then the transformation matrix used is 


p q 



which is orthonormal. For simpler notation, C - cos G, and 
S •- sin 6 has been used. 
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r 

The value of 6 Is calculated as shown In the following ex- 
ample. A system of order 5 with pivot location p « 2, q 4 Is 
used, but the results may be applied to any size system with any 
pivot location p, q. 


■ < 


lAl,. 

1 

k 














r 1 




"u 

“ip 

®13 

“iq 

*15 


1 




“pi 

“pp 

“p3 

“pq 

“p3 


C -b 




^31 

“3p 

^33 

“3q 

*35 


1 




“qi 

a 

qp 

“q3 

“qq 

“q5' 


s c 




/31 

“5p 

*53 

“3q 

*53_ 


1 








a, <■ 

Ip 

: 4 

• V 

*13 


-AjpS ♦ *15 

. + a ,S 
qi 


a 

PP 

+ 2a CS > a S*- 

pq qq 

a -C 4- a ,S 
p3 q3 

l“qq 


^31 




a. C 
3p 

+ 

S 

3q 

*33 


-*3p® * ‘‘3.,'^ “35 


a ,S + a ,C (a - a ICS + a (C*’ - S^) -a .S + a ,C 

pi <ji \ qq pp/ pq pj 

^51 *3p^ ^5q^ ®33 


a S-^ - 2a CS + a 
pp pq qq 

-a, S + a, C 
3p 3q 


-a ,S + a .C 
p5 q5 




( 6 ) 


where use is made of the symmetry property a 


pq 


seen, only the p, q rows and columns of [A],, 
calculating [A]j^. 


= a . As can be 

qp 

are altered in 


Setting the pivot element p, q of [A]j^ to zero gives 


a (cos^ 0 - sin^ 6) * /a - a \ sin 6 cos 6 

pq \ PP qqj 


from which 

2a 

tan 0 = , ~ (7) 

la - a \ ± Jla - a \^ + 4a^ 

I PP qqj t1 pp qqj pq 

The sign used with the radical is the same a& the sign of 
(a - a \. This is done to avoid small differences of large 

\ pp qqj 

numbers. 
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Using trigonometric identities gives 

cos 0 = i/Vr + tan^ 0 (8) 

and 

sin 0 « tan u cos 0 (9) 

as the values to use in Equation (6) to set a =0. 

pq 

Tt should be obvious that even though the pivot element has 

been reduced to zero by one rotation, it may take on a non-zero 

value in some later rotation. A rotation gives a step tov’ rd the 

diagonal form because the sum of the squares of the off-diagonal 

elements of [A], is less than that of [A], , by the amount 2a** . 

k k-1 ^ pq 

The importance of the threshold version is re-emphasized in that 

small values (relatively) of the off-diagonal elements of [A] are 

not used for pivot elements. Proof of convergence is in the 

literature (References 2 and 3) and will not be given here. 


As mentioned before, the continued product [T]^ [T]^ [1]^^ 

converges to the eigenvectors [4>]. Thus, would be calculated 






m 

"^ll 

"^12 

'^13 

^4 


"^21 

^*22 

♦ 23 

'*'24 


'^31 

1 

‘^32 

4-33 

'*'34 


1 

!^i 

"^42 

'^43 

^4 


Si 

^2 

♦53 

^4 

m 

”♦11 

<P 

12^-^ 

'*'14^ 

i 

*21 


22C + 

♦24S 


'*’31 


32C + 

«34S 




42^-^ 

^44^ 


Si 


52^-^ 

♦54S 


♦ 15 “ 


1 

^25 


c -s 

♦35 


1 

^5 


s c 

^5_ 


1 


13 


+ ^4<^ 

•hs 

23 

••«22S 

+ ^24^ 

‘*’25 

33 

-♦32=" 

+ 

♦35 

*43 

-*42S 

+ */,4C 

'■45 

'53 

••^52^ 

+ ♦saC 

^5 


( 10 ) 
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for our previous example of a system of order 5 and pivot location 

p = 2, q = 4. As can be seen, only the p, q columns of [‘t] are 

K*! 

altered in calculating 

The resulting eigenvalues and eigenvectors can be checked as 

T 

follows. Post multiplying Equation (1) by [tl*] gives 

[A1 = [<f] t X ] (11) 

because the eigenvectors are orthonormal. Therefore, to check the 

T 

results, perform the triple matrix product t ^ ] l^] and com- 
pare with the original [A] . This operation is not performed in the 
subroutine. 




are 
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and 

['Ml = 

/2T3 

1/i^ 

0. 



-Ih^ 


0. 



0. 

0. 

1. 


Second pivot ; 


p = 2, q = 3 


“p., ■ “23 ■ 


i == = .5 

PP 22 


qq 33 


From Equations (7), (8), and (9) 
tan e = I //3 
cos G = 
sin 0 = .b 

From Equations (6) cind (10), the results of the second pivot 

are 


[M 2 = 

0 

0. 



0. 

0. 



(sym) 

2. 


[$l2 = 


.5 

i_7' 

~ 


1 

“ 

1 

/2 

1 

" /6 


0. 

.5 

Jz 

2 


The diagonal elements of the [A ] 2 are the eigenvalues *^nd the 
volumns of corresponding eigenvectors. 
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Subroutine EIGNIA calculates the eigenvalues (proper values , 
characteristic roots, latent roots) and eigenvectors (proper 
vectors, characterise i.c vectors, latent vectors) of a real sym- 
metric matrix using method of C. G- J. Jacobi. This subroutine 
is a modification of Subroutine EIGNl and uses a convergence 
tolerance on X instead of a final off-diagonal value to terminate 
the solution. The eigenvalue problem can be expressed as 

[A] [4-1 = M W [11 

where 

[A] = a real symmetric matrix of order N; 

[($] = a matrix whose columns are the eigenvectors of Eq [1]; 

fXj = a diagonal matrix whose elements are the corresponding 
eigenvalues . 

This method diagonalizes matrix [A] by successive ^'lane rotations 
and leads simultaneously to all eigenvalues and eigenvectors. 
Theoretically, an infinite number of plane rotations are necessary 
CO produce the diagonal form with all off-diagonal elements equal 
to zero. Practically, the number of plane rotations is limited to 
a finite number by stopping the process when [X^(k) - X^(k-1]/ 

X^(k), i ® 1, N is less than a prescribed value (k refers to 

plane rotation number). This prescribed value (denoted as CTVAL 
in the subroutine) may be selected and input by the analyst or 
at the analyst's option determined in the subroutine. 

A threshold version is used» that is, each off-diagonal element 
of (A] is compared in regular sequence with a threshold value 
and a plane rotation performed only if its magnitude exceeds the 
threshold value. Tnis is done to give faster convergence to the 
diagonal form. Because [A] is symmetric, only the upper half 
needs be examined. The threshold value is lowered whenever there 
are no remaining off-diagonal elements with magnitudes larger than 
the threshold. 

The initial threshold value in the subroutine is obtained by divid- 
ing the maximum off-diagonal element of [A] by 10. Subsequent 
reductions in tbp threshold value are accomplished by dividing the 
threshold value by 10 also. 

The threshold version is the fastest version of Jacobi's method 
(Ref 2). 

An important feature of the Jacobi method is that equal eigenvalues 
and the corresponding eigenvectors are easily calculated with no 
change in the algorithm. 
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DESCRIPTION OF TECHNIQUE 

Assume that an order matrix [T] exists that in orthonormal » 
that is, (IT]^ IT] = [iD such that 


(T]"^ [A] IT] = tDj. 12] 

where fDj is a diagonal matrix and [A] is the original matrix. 
Premultiply Eq [2] by [T] to give 

[A] [X] = [Tl TDJ. 13] 

Comparing EQ [1] and [3], it is seen that if matrix [T] can be 
found which transforms a real symmetric matrix [A] into a diagonal 

matrix roj* then the i^^ diagonal element of fDj is the i^^ eigen- 
value of [A], and the i^^ column of [T] is the i^^ eigenvector of 
lA]. 


The orthonormal matrix [T] is built in stepwise fashion using 
elementary orthonormal transformations [T]^ to annihilate* in 

turn* selected off-diagonal elements of [A]; that is, starting 
with the given matrix [A], operate as follows* 


lA] 


= [A]q 


IT]\ [A]q IT]^ = [A]^ 

IT]^ [A]^ [T]2 = [A]2 


IT] 


T 

k 


[A] 


k-1 


lT]k 


l^lk 


SO that 


[T]]^ . . . [T]^ IT]^ lA] [T]^ (T]2 . . . [T]j^ = [A]^^. [4] 

Comparing Eq [4] with Eq [2], [T]j^ must be chosen so that the 
triple matrix product result, (Ajj^ converges to fDj = fAj and the 
continued product [T]^ [T]^ ••• [T]j^ converges to [T] = The 

most important property of is that each triple matrix product 

(rotation) with [A]j^ ^ in Eq [A] causes a particular off-diagonal 
element in [A]^ ^ to vanish. As mentioned before, an off-diagonal 
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element of [A]j^ ^ is picked for elimination if It exceeds a thres- 
hold value and is referred to as the pivot element. If the p, q 
element of [A]. . meets the criteria, then the transformation matrix 
used is 



which is orthonormal. For simpler notation, C= cos 6, and 
S5 sin 6 has been used. 


The value of 6 is calculated as shown in the following example. 
A system of order 5 with pivot location p = 2, q = 4 is used, 
but the results may be applied to any size system with any pivot 
location p, q. 
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where use is made of the symmetry property As can 

be seen, only the p, q rows and columns of [A], . are altered in 
calculating 

Setting the pivot element p, q of lAjj^ to aero gives 
2 2 

a (cos 9 - sin 6) = (a - a ) sin 9 cos 0 

pq IP qq 


from which 


2a 



The sign used with the radical is the same as the sign of 
(a - a ). This is done to avoid small differences of large 

pp qq 

numbers. 


Using trigonometric identities gives 

cos 6 “ 1/.0 + tan^ 6 
and 

sin 6 = tan 9 cos 0 

as the values to use in Eq [6] to set a = 0. 


[ 8 ] 

19] 


It should be obvious that even though the pivot element has been 
reduced to zero by one rotation^ it may take on a nonzero value 
in some later rotation* A rotation gives a step toward the 
diagonal form because the sum of the squares of the off-diagonal 
elements of less than that of by the amount 

2a^ • The importance of the threshold version is reemphasized 
P<1 

in that small values (relatively) of the off-diagonal elements of 
[A} are not used for pivot elements. Proof of convergence is in 
the literature (References 2 and 3) and will not be given here. 


As mentioned before, the continued product (Tjj^ [TJ 2 ••• [Tjj^ 
converges to the eigenvectors Thus, would be calculated 





m 


♦12 

♦l3 

♦l4 

♦is 


I 




♦22 

♦23 

♦24 

♦25 


c -s 



♦31 

♦32 

♦33 

♦34 

♦35 


1 



♦41 

♦42 

♦43 


♦45 


s c 



♦si 

♦S2 

953 

♦54 

♦ss_ 


1 

■ 

m 

"♦ll 


^IjC + 

♦ 14S 

♦l3 


♦is" 


♦21 

»22C + »24S 

♦23 

-«22S + 

♦25 


♦31 


932C + ♦34S 

♦33 

-♦32" + *34'= 

«3S 


♦41 


♦42^^ * ^44^ 

♦43 

■♦42® + ^44*^ 

♦45 


♦si 


♦ 52C + 

♦54^ 

♦53 

-♦52® *54^ 

♦55 
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tioi 


for our previous example of a system of order 5 and pivot location 
p ■ 2, q = As can be seen, only the p, q columns of [$] are 

altered in calculating 

'I'he resulting eigenvalues and eigenvectors can be checked as 
follows. Postmultiplying Eq [1] by [4>]^ gives 

[A] » t4>] fAj 1$]^ 

because the eigenvectors are orthonorraal. Therefore, to check the 

results, perform the triple matrix product [4>] fXj [<t]^ and com- 
pare with the original [A]. This operation is not performed in 
the subrou*-ine. 


Example 

Consider [A] e 


1.5 -1/*^ -0.5 

1.0 - 1//2 

(sym) 1.5 
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First pivot: p » 1, q = 2 ’ 


*pq * ^12 “ 

%P “ hi “ 
*qq “ ^22 ” 


From Eq (71. [81, and [91 

tan 6 « -Hy/i 
cos 9 « »T/T 
sin e = -l/i/J 

From Eq [61 and [101, the results of the first pivot are 


(Alj^ « 2. 0. 


0.5 -✓372 


(sym) 


^ 1//3 0. 

- 1//3 ^In 0 . 


Second pivot: p » 2, 1=3 

‘n " “23 ■ 

*pp ■ “22 ■ 

“qq ■ “33 ■ 

From Eq [7J, [8], and [9] 

tan 0 = 1//Z 


cos 0 = ✓3/2 


sin 9 = 0.5 
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From Eq [6] and [10] , the results of the second pivot are 


IAI2 


2 . 


0 . 


0 . 


0 . 0 . 


and 

= 


jj[sym) 




0.5 

_ 1 

1 

Jl 

_ 

' /I 

Jl 

' ^ 

0 . 

0.5 

/3 

2 


The diagonal elements of the [A]^ are the eigenvalues and the 
columns of [$]2 ^te the corresponding eigenvectors. 
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Subroutine FRl calculates the frequency response {X(w)} from 
the following differential equation expressed in the frequency 
domain: 


(-oj 2(A] + iw[B] + [C]) {X(w)> = {D]{F(6 u)}. (1) 

Matrices [A], [B], [C], and [D] are input directly to this sub- 
routine along with a set of values of w. For a structural prob- 
lem, [A] is the mass matrix, [B] the damping matrix, [C] the 
stiffness matrix, and [D] is the transpose of the vibration mode 
shapes or a unity matrix depending on whether the equations are a 
modal or discrete representation of the structure. The force 
{F(ui)}, assumed real, is calculated internally in the subroutine 
using linear interpolation ^^rith [TABW] and (TABFl which are both 
input to the subroutine. As an illustration of the use of [TABW] 
and [TABF] consider the following example: 


Fl 


2 . 0 # 




1 ■■■ ■ .. » II ■■ ■ ^ U) 

0 5 10 15 20 


1.0 


F2 


.5 



°0 5 10 15 20 




The table (A table is defined in this report as a matrix that 
may have incomplete column data in some rows.) giving the independ- 
ent variable coordinate w is 


[TABW] = 


0 . 


20 . 


5. 10. 15. . 




The table giving the corresponding coordinates of the dependent 
variable F is 


[TABF] = 


2 . 


2 . 


.5 1. 


.5 
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The following values for {F(i»/)} will be obtained by linear inter- 
polation at w = 7. 


{F(w = 7)} 


2 . 


.7 


DESCRIPTION OF TECHNIQUE 

First, consider the inversion of a complex matrix. Let ([G] 

+ i[H]) be the original matrix expressed as a complex matrix in 
terms of two real matrices [G] and_[H] . The product of a matrix 
with its inverse (assumed to be ([G] + i[H])) is, by definition, 
unity; thus 

(IG] + i[H])([G] + i[H]) = [G][G] - -i- i(tH][G} + [GHH]) = (I] + i[0] (2) 


or, equating real and imaginary parts, 

[G] [G] - [H][H] = [I] (2a) 

[H] [G] + tG](H] = [0]. (2b) 

Two approaches now exist for the solution. From Equation (2b), 

either [H] = ~[G][H][G]~^ (assumes [G] is nonsingular), (3a) 

or [G] = -[H][G][H]“^ (assumes [Hj is nonsingular). (3b) 

Inspection of Equation (1) shows chat 

[G] = -[.C] - (o^[A] 
and [H] = w[B]. 


Because many values of a> can make [G] singular (and also ill-con- 
ditioned such that the inverse is bad), it was decided to use 
Equation (3b) which then requires the damping matrix [B] to be 
nonsingular. Substituting Equation (3b) into (2a) yields 

[H] = -([G] [H]-1[G] + [H])";! (3c) 

To eliminate having to do the triple matrix product [G][H]“'[G] for 
every value of o), Equation (1) is modified so that [H] = oj[lj. This 
is done by premultiplying by [B]~^ so that 


(-oj^tB]-llA] + im[l] + [Bj-ltC]){X(w)} = [B]-l[D)lF(oj)] . 


(la) 
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Solving for {X(w)} gives 

{\eal^‘"^} “ tG][B]’l[D]{F(aj)} (4a) 

and * I“HB]-1 [D]{F(u))} (4b) 

where 

[G] “ “ ^ fB] IG] (4c) 

[H] = -oj([G] 2 + 0)2[I])-1 (4d) 

and [G] * tB]-l[C] - u)2[B]-l[A]. (4e) 
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Subroutine FRAEl calculates frequency response additional 
equations. That is, 

{Z(o))} == (A]{X(u))} (1) 

where {X(u))} is the complex response previously calculated and 
written on tape in frequency response subroutine FRl. That tape 
is read in this subroutine FRAEl to calculate {Z(u))}. 

If {X(o))} is the response calculated from a modal representa- 
tion of the structure, then [A] would be the mode shapes so that 
Equation (1) would then give {Z(w)} as the response of the dis- 
crete system. {Z(u>)} will also be com[;iex. {Z(w)} is printed for 
every w value of {X(to)} from tape NXXAPE and also written on tape 
HZTjsFE for any subsequent use. 
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Subroutine INTAPE initializes a tape (disk is preferred, see 
writeup of Subroutine WTAPE) for the FORMA tape system by writing 
EOT (end of tape) at the beginning of the tape (disk). AaI FORMA 
tape subroutines recognize this EOT as being the end of written 
data* Each **rew** tape (disk'^ must be initialized with this Sub- 
routine INTAPE to make the tape (disk) compatible with the other 
FORMA tape subroutines (LTAPE, RTAPE, WTAPE, and UPDATE). 

A "new'* tape (disk) is defined as a tape (disk) for which it 
is desired to start writing matrix data at the front of the tape 
(disk)* Thus, a "new" tape (disk) could be one with obsolete 
FORMA matrix data on **.t as well as one that has never been written 
on by the ^^ORMA system* 

As an example, pertinent statements from a program containing 
INTAPE could be: 

DATA NIT, NOT/5, 6/ 

1001 FORMAT (12A6) 

NRTAPE = 10 

READ (NIT, 1001) IFINIT, TAPEID 

IF (IFINIT .EQ. 6HINITIL) CALL INTAPE (NRTAPE, TAPEID) 


The input data (starting in card column 1) to this example 
program would be: 

either INITILTXXXX, if the tape is to be initialized. 

(TXXXX represents the particular tape number 
used, e.g«, T1234); 

or NOINIT, if the tape is not to be initialized. 

The tape Identification is not needed. 
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Subroutine INVl inverts a matrix by the bordering method (Ref 
1) . In matrix notation, the inversion is represented by 

where 

[A] is the matrix to be inverted, 

[Z] is the result, and 
N is the size of the matrices. 

Matrix [A] may be nonsymmetric. The determinant ratios of [A] are 
also calculated and printed. The determinant of [A] is the prcduct 
of these determinant ratios. 

An l?.version check, [Z] [A], is calculated (The product should 
be a unity matrix.) and a summary of these results is printed. 

REFERENCE 

1. Bodewig, E.: Matrix CaiouluSy North Holland Publishing Com- 

pany, Amsterdam, 1959. 
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Subroutine INV2 inverts a matrix by the rank annihilation 
method (Ref 1). In ?.atrix notation, the inversion is represented 
by 


where 


[A] is the matrix to be inverted, 

[Z] is the result, and 
R is the size of the matrices. 

Matrix [A] may b'* nonsymmetric. 

An inversic • ,.k, [Z] lA], is calculated (The product should 

be a unity matrix.) and a summary of these results is printed. 

REFERENCE 

1. Ralston, A., and Wolf, H. S.: Mathematical Methods for Digital 

Computers, John Wiley & Sons, 1966. 



Subroutine INV3 Inverts a symmetric matrix using triangular 
decomposition and triangular inversion « In matrix notation > the 
inversion is represented by 


where 


[A] is the symmetric matrix to be inverted > 

[Z] is the result (also symmetri :) » and 
N is the size of the matrices. 

In addition to being symmetric, matrix [A] should be positive 
definite. Because symmetry is assumed, only the upper half of 
[a] will be used to calculate [Z] which is also symmetric. The 
determinant ratios of [A] are also calculated and printed. The 
determinant of [A] is the product of these determinant ratios. 

An inversion check, [Z] [A], is calculatea (xhe product should 
be a unity matrix.) and a summary of these results is printed. 

Only the upper half of [A] is used to calculate [Z] , but all of 
[A] (and [Z]) is used in the inversion check. Thus» errors in 
the inversion check can result from a nonsymmetric [A] as well as 
from a poor inversion. 

DESCRIPTION OF TECHNIQUE 

The first operation in the inversion of [A] is to decompose 
[A] (reference Subroutine DCOMl) into triangular factors, that 
is. 


[A] = [U]^ [U] 


where [U] is an upper triangular matrix. Then 


[Z] - 


The Inversion of the upper 
by Subroutine INV4. 


[A]-l 


([Ul" 

lul)-‘ 



tul'i 

(lurf. 

triangular matrix 


[U], is calculated 
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Subroutine INV 4 inverts an upper triangular matrix. In matrix 
notation the inversion is represented by 

where 

[Al is the upper triangular matrix to be inverted, 

[Z] is the result (also upper triangular) , and 
N is the size of the matrices* 

An inversion check is not calculated* 

DESCRIPTION OF TECHNIQUE 

That [Z] is upper triangular if [A] is upper triangular is 
easily demonstrated by the following example* From [A] [Z] == [I], 


an 

ai2 

ai3 


211 

212 

213 



0 

0” 

0 

S22 

32 3 


Z21 

222 

22 3 

= 

0 

1 

0 

0 

0 

333 


Z31 

232 

233 


0 

0 

1 


Preforming the matrix product and equating like elements gives 

( 3,1 element) a33 Z31 =0 

Because a33 # 0 , Z3X = 0 * 


(3,2 element) 


^33 

Because a33 


Z32 = 0 

0 *•* Z32 = 0. 


(2,1 element) a22 ^21 ^23 ^31 “ ^ 

Because Z31 = 0 and a22 ^ 0 , *•. Z21 “ 0 . 
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The algorithm used to calculate [Z] Is obtained by consider- 
ing the equation [A] [Z] » [1] as shown below for size N « 4. 




Equating elements of the matrices on the left and right side of 
the equal sign gives 

*11 “ V*li 

212 = “ 2X2 Z 22 / 21 I 

= - zxi ai2 Z 22 

223 = “ 223 Z 33/222 
= - Z22 323 233 

Z34 = “ 234 Z44/a33 
“ “ 233 334 Z44 

213 * ~ (212 Z23 + 213 233)/aii 

= - (zii ai 3 + 2 x 2 323) Z33 

224 “ ~ (223 234 + 324 Z44)/a22 

= - (z 22 324 + Z 23 334) Z44 

Zx4 “ ~ (ax2 224 + ax3 Z34 + 3X4 244)/axi 

= - (zxi 3x4 + ZX2 324 2x3 334 ) Z 44 
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From the above equations, the algorithm for the (i,j) element of 
[Z] for size N is given in general as 


z 


ii 



(i = 1,N) 


j-1 

E 

k=i 


' '"JJ La "ik "kj 


/i = 1, N - 1 

U > i 



LTAPE 


Subroutine LTAPE lists the matrix headings (see Subroutine WTAPE 
and YWTAPE writeups) written on a FORMA tape (or disk) . These 
matrix headings were written by Subroutines l^TIAPE and YWTAPE and 
consist of : 

NO. » Matrix number on tape, 

RUN NO. = Run number of problem when matrix was written on tape, 

NAME = Matrix name; 

NROWS = Number of rows of matrix; 

NCOLS = Number of columns of matrix; 

DATE = Date when matrix was written on cape; 

NNZ = Number of nonzeros (used only in sparse FORMA when 

only nonzeros are used) ; 

PARTITION = Partition number of sparse matrix. 
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Subroutine MASSl takes (on option) distributed nass> distrib- 
uted rotary inertia, and concentrated mass items of a beam and 
replaces them with a mass matrix. The elements of the mass matrix 
are representative inertial values of the beam at selected points 
on the beam. These elements are calculated by assuming a linear 
function between pairs of the selected points to describe the 
beam's lateral elastic deflection. 

The x-stations of the selected points (panel points) are 
given in {PP}. These x-stations must be in increasing order. 

The distributed mass, m(x) , is assumed to be piecewise linear 
and is represented by straight line segments as shown in Figure 1 . 
The x-stations of the end points for the line segments giving the 
distributed mass are independent of the panel point x-stations. 
However, the distributed mass must be within the panel point lim- 
its. The line segments representing the distributed mass may or 
may not be joined and may overlap. The distributed mass is de- 
fined in [DMASS]. Each row of [DMASS] represents one nonvertical 
line segment. The form of each row of [DMASS] is [xj X2 m^ m2] 
where x^ , mj give the first end point and X2, m2 give the second 
end point of a line segment. 


m(x) 



The distrlDuted rotary inertia, r(x), is also assumed to be 
represented by straight line segments. All statements used above 
for distributed mass are applicable to distributed rotary inertia. 
The distributed rotary inertia is defined in [DRIM] . The form of 
each row of [DRIN] is [xj X2 rj r2]. The end point stations X] 
and X2 for the rotary inertia line segments are independent from 
the end point stations and for the mass line segments. 
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The concentrated mass items are defined In [CONC] • Each row 
of [CONC] represents one concentrated mass item and contains: 


1 ) 

2 ) 

3) 

4) 


X 9 the station at which the xtem is attached to the 
a 

beam; 

M , the mass of the item; 
c 


X , the center of gravity of 
eg’ 

eg 

I the moment of inertia of 
c 

center of gravity* 


the item; and 

the item about its own 


These four elements are given in the form x M x . The 

L a c eg c J 

attachment station, x , of an item to the beam must be within the 
panel point limits. 


The calculated representative inertial values at the selected 
panel points are placed in a mass matrix given by [Z] which is 
symmetric and tridiagonal* That is. 


^ 1.1 ^ 1,2 

*2,1 *2,2 *2,3 


UJ 


n»n 


*3,2 *3,3 

^ 4,3 


'^ 3.4 

""4,4 *4,5 


L *n-l,n-2 *n-l,n-l ^n-l,n 

z , z 

n,n-*l n,n 

where z, , = z, . and n is the number of panel points. The gen- 
i > J 1 » i 

eralized coordinates associated with [Z] are lateral translation 
at each of the panel points* 
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The sign convention used in this paper to obtain the mass 
matrix [Z] is shown in Figure 2. 




0 


" I 

Figure 2 Sign Convention for This Paper 


The beam mass matrix obtained in this paper is applicable for 
either the pitch, yaw, roll, or longitudinal plane with a change 
of variables. For the axis system shown in Figure 3, the follow- 
ing variables would be used: 




Axis System of Figure 3 

This 
Paper 
(6, x) 

Pitch 

Yaw 

Roll 

Longitudinal 

(z, x) 

(y. x) 

X, 0 

X 

(x) 

X 

X 

X 

X 

X 

6 

6 

6 

e 

6 


z 

y 

X 

X 

m(x) 

m(x) 

ia(x) 

i,(x) 

m(x) 

r(x) 

r (x) 

y 

r^(x) 

NA 

NA 

M 

M 

M 

1° 

M 

c 

c 

c 

c,x 

c 

iCg 

c 

icg 

c,y 

iCg 

c,z 

NA 

NA 

1 NA » Not Applicable 




A z 





X 


u 

X 


Figure 3 Beam Axis System (Right Hand) 
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The beam's total mass, center of gravity, and moment of iner-* 
tia about the beam center of grc^.ity are also calculated and 
printed. 


DESCRIPTION OF TECHNIQUE 

The replacement of distributed mass, distributed rotary iner~ 
tia, and concentrated mass items of a beam by a mass matrix is 
obtained using a kinetic energy approach as follows. 


For small deflections normal to the longitudinal body axis 
(x) of the beam shox^n in Figure 2, the kinetic energy is defined 
by 


where 


T = 


E 

h f [m(x) 6^(x,t) + r(x) 6^(x,t)] dx 


( 1 ) 


m(x) is the distributed mass, 

r(x) is the distributed rotary inertia, 

M is the mass of a concentrated item, 
c 

CR 

I ^ is the rotary inertia of a concentrated item about 
^ its own center of gravity, 

6 is the time rate of change of elastic deflection 
normal to the body x-axis, referred to as lateral 
velocity in this paper, 

6 Is the angular velocity about an axis y^i^pendlcular 
to the paper, 

A is the velocity of a concentrated Lwem and will be 
explained later, 

X is the undeformed longitudinal axis of the beam, 
t is time, 

Xg is the starting x-station of the beam, and 

Xg is the ending x~station of the beam. 

The finite summation is over the iiu«.ber of concentrated mass 
items. 
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The technique used here to integrate Equation (1) is es- 
sentially a modification of the Rayleigh-Ritz method. In the 
basic Rayleigh-Ritz method the elastic deflection of a beam is 
represented by a summation of assumed displacement functions each 
weighted by a generalized coordinate representing the contribu- 
tion of each of the assumed functions. Each displacement func- 
tion is assumed over the entire beam len^^’ and considerable 
skill is required to choose these functions. The modification 
used here is to assume a simple displacement function between 
consecutive panel points and to use the panel point displacements 
as the generalized coordinates. A linear displacement function 
will be assumed here. This is illustrated in Figure 4 where 

and consecutive panel point stations. The region between 

panel points k and k+1 is referred to as bay k. The lateral dis- 
placement between panel points k and k+1 is given by 

6(K) = \ + (« - \) (\+l - \)/(*k+l - »k)- 

Similarly, the lateral velocity is given by 

6(x.t) - ( X - - 6 Jt)] I - x^j. (2) 

The angular velocity is obtained as the geometric derivative of 
the lateral velocity. That is. 


9(x,t) 


36(x.t) 

9x 


Pk+l^*^^ ~ I(\+1 ' 


( 3 ) 


This equation for the angular velocity is slightly in error be- 
cause a rotation due to shear deformation is included. A knowl- 
edge of the stiffness distribution would be necessary to deter- 
mine the relative contribution of bending rotation and shear 
deformation to the angular velocity. This error will be small 
for beam type structures and will be ignored# 
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The kinetic energy of distrijLted mass, distributed rotary 
inertia, and concentrated mass items will be considered separately* 
The distributed mass, m(x) , is considered first. The geometry 
for a line segment of distributed mass is shown below in Figure 

5* 



The equation for a straight line segment as shown in Figure 

5 is 


m(7i) - mi + (x “ xi) (m 2 - mi)/(x 2 “ xi), (4) 

Substituting Equations (2) and (4) into (1) gives the lecic 
energy of the distributed mass represented lint, *>egment 

i in bay k as 
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The subscripts p and q have been introduced to handle the possi- 
bility of a line segment extending past the bay limits. Thus, x 


is the greater of or 


and X 

q 


is the lesser of xo or x, . . • 
^ k+1 


Similarly, m^ is either m^ or m^^ and m^ is either m2 or 

The integration is continued for the line segment in adjacent 
bays, if necessary, until the entire line segment has been used. 
Performing the integration of Equation ( 5 ) yields 




6,(0 






k,k 

(sym) 


^k,k+l 

^k+l,k+J V 


\(t) 1 




where 


Fl - 2F2 + F3 

\,k+l = ^2 - F3 

\+l,k+l 

^ ■ %i"p + ”,1/2 

F2 * L l/m + m \ /h + H \ + m H + m H I /6 
p |( p q) \ p q/ p p q ql/ 

Fq * L 


L - X - X 

p q p 


/m + m V m + H + 2 (m + m H^\\ 

Ip q) \ p 0/ \ p p q q/j/ 

)A 


X - X, 

p k. 


“p- ( 

■ C*, ■ \A' 

\ ■ \+x • 


( 6 ) 


(6a) 

(6b) 

(6c) 

(6d) 

( 6 »- 

( 6 £) 

(6g) 

(6h) 

(6i) 

(6j) 


The kernel matrix in the triple matrix product of Equation 6) 
is the mass matrix that replaces the distributed mass repr -rented 
by one line segment i in bay k. 
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The kinetic energy of the distributed rotary inertia, r(x), 
is considered next. The geometiry for one rotary inertia linr 
segment is similar to that given in Figure 5 for distributed mass. 
Similarly y 


r(x) = ri 4- (x - xi) (r -2 - ri)/(x 2 “ xj). (7) 

Substituting Equations (3) and (7) : to (1) gives the kinetic 
energy ^^f the distributed rotary intitia represented by one line 
segment i in bay k as 



P 


[- - 5^(0) /(vi - \)] <8) 

The subscripts p and q have similar meaning as has been previously 
discussed for distributed mass. 


Performing the integration of Equation (8) vields the sair.^ 
equation as was obtained previously for distributed mass. It is 
repeated here as 




'k,k+l 


k+l,k+lJL k+1 


6k 't) 1 


(t) 


(9) 


Now, however. 


5 , 1, = z, . 1 , . 1 = L /r + r \I 2 L^ 

k,k k+1, k+1 p\ p q// k 


(9a) 


and 


\,k+l " "\,k* 


As before. 


L =• 
P 


X 

q 



(9b) 


(9c) 


and 




^k ■' ^k+1 


(9d) 
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The kernel matrix of Equation (9) is then the mass matrix that 
replaces the distributed rotary inertia represented by one line 
segment i in bay k. 

The CO . 4 rntrated mass items are considered last. These items 
X- ;uire investigation because the item may be connected to 

the beam at an attach point a by an arm (assumed rigid) as shown 
in Figure 6 * The square of the velocity of concentrated item c 
is 

4^(0 - [4^(0 - - xj + [■>%<'>]". «“> 



From Equation Cl) » using the velocity expression of Equation 
( 10 ) » the kinetic energy for a concentrated item c becomes: 

[«c - »a) °a‘'> * “I'c ‘ ’‘a)'“ 

( 11 ) 


• ♦ 

where use has been made of the relationship that 6 (t) = 6 (t)» 

c a 

In further work, the offset D will be assumed to be zero. The 
reason for this is that any mass offset from the beam x-axis means 
that the longitudinal principal axis of the beam does not coincide 
with the beam x-axis* A coupled transverse-longitudinal analysis 
would be required to accurately describe the kinetic energy of the 
beam for D 5 ^ 0. However, if the analyst wishes to include the ef- 
fect of an offset D, inspection of Equation (11) shows that 

4* M could be used. With the above conditions, the kinetic 

c c 

energy for a concentrated item can be expressed in matrix notation as 
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( 12 ) 


The lateral and rotational velocities at the attach point a are 
given in terms of the velocities at adjacent panel points k and 
k+1 from Equations (2) and C3) > respectively. Using these equa- 
tions the kinetic energy of one concentrated item c in bay k is 
finally given as 




(t) 



^k,k+l ^ 

'k+l,k+lj 


,(t) 

k+1^*^^ 


where 


(13) 


^k,k 


= F, 


2F" + Fq 


(13a) 


\,k+l 


= F2 


F 3 , and 


(13b) 


\+l,k+l 


= F3, 


(13c) 


as were previously obtained for 'istributed mass. Now, however, 

Fi = M , (13d) 

C13e) 
(13f) 
(13g) 


‘■k " - *k- 

The terms in Eqiations (13a) and (13b) could be combined to give 
a simpler expression for and k+i* used, however, 

to minimize coding instructions in the subroutine. The attach 

point coordinate, x , does not appear in any of the final equations 

a 


Fo = M H 
^ c c 


Fq = M 


h 2 + 

c c / b 
"c ' (*c - \)jh- 


and 


As before. 
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for etc. It is used only to determine the bay in which the 

attach point is located. For an attach point coinciding with a 
panel point, the inertial properties could be correctly calculated 
using the coordinates of the panel point bay in front of the at- 
tach point or the panel point bay in back of the attach point - 
The results would be different, however, since the angular veloc- 
ity is discontinuous at a panel point. To minimize numerical 
roundoff error, the attach point is assumed at the beginning of 
the bay. 

The mass matrix for the entire beam is obtained by evaluating 
Equations (6) for each line segment of distributed mass. Equations 
(9) for each line segment of distributed rotary inertia, and 
Equations [13] for each concentrated mass item. Every resulting 
2x2 bay mass matrix is added to previous 2x2 bay mass matrices at 
like panel points to form the mass matrix for the entire beam. 

The total mass properties of the beam are calculated by the 
following triple matrix product. 



where 

{1} is a column of ones, 

{PP} is a column of the panel point x-stations, 

[Z] is the mass matrix for the beam, 
is the beam mass, 

o 

is the first moment of the beam about x = 0, and 

is the moment of inertia of the beam about x » 0. 

From this data the center of gravity of the beam is calculated 
from 



and the moment of inertia of the beam about x is calculated from 

eg 
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SPECIAL CASES 


Two cases are worthy of mention. In the first case* a Uncai'li 
varying segment of distributed mass extending to the bay limits 
|Xj^ and has a mass matrix (from Equations (.6) of; 


^k-i-1 \ 


“k+i "H. * "k+i 


(sym) 


"he . 


In the second case, a unifoim |i.e., nij^ = ~ segment of 

distributed mass extending to the bay limits /x and x, .. j has a 
mass matrix of : ^ 




MASS 2 — 1/16 


Subroutine MASS 2 takes (on option) dxstribuiied mass, distrib- 
uted rotary inertia^ and concentrated mass items of a beam and 
replaces them with a mass matrix- The elements of the mass matrix 
are representative inertial values of the beam at selected points 
on the beam. These elements are calculated by assuming a cubic 
function between pairs of the selected points to describe the 
beam's lateral elastic deflection. The x-stations of the selected 
points (panel points) are given in IPP). These x-stations must 
be In increasing order. 

The distributed mass, m(x) , is assumed to be piecewise linear 
and is represented by straight line segments as shown in Figure 1 . 
The x-stations of the end points for the line segments giving the 
distributed mass are independent of the panel point x-stations - 
However, the distributed mass must be within the panel point lim- 
its. The line segments representing the distributed mass may or 
may not be joined and may overlap. The distributed mass is de- 
fined in [DMASSj. Each row of [DMASS] represents one nonvertical 
line segment. The form of each row of [DMASS] is [xj X2 mj m2] 
where , mi give the first end point and X2, m2 give the second 
end point of a line segment. 


m(x) 



Figure 1 Distributed Mass 


The distributed rotary inertia, r(x), is also assumed to be 
represented by straight line segments. All statements used above 
for distributed mass are applicable to distributed rotary inertia. 
The distributed rotary inertia is defined in [DRIN] . The form of 
each row of [DRIN] is [xi X2 ri r2]. The end point stations Xj 
and X2 for the rotary inertia line segments are independent from 
the end point stations xi and X2 for the mass line segments. 
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The concentrated mass items are defined in [CONC] • Each row 
of [CONC] represents one concentrated mass item and contains: 

1) X . the station at which the item is attached to the 

a 

beam; 

2) M , the mass of the item; 

c 

3) center of gravity of the item; 

4) the moment of inertia of the item about its 
c 

own center of gravity. 

These four elements are given in the form x M x . The 

a c eg c j 

attachment station, x^, of an item to the beam must be within the 
panel point limits. 


The calculated representative inertial values at the selected 
panel points are placed in a mass matrix given by [Z]. The form 
of [Z] is 


I^l2nx2n 

where n is the number of panel points. Matrix [Z] is S 3 rmmetric, 
i.e., = 2 ^^. The partition form of [Z] results from the two 

generalized coordinates (lateral tj.anslation, 6, and rotation, 6) 
at each panel point arranged with all translation coordinates 
first followed by all rotation coordinates. Each partition of 
[Zj is square and tri-diagonal, for example, 

1.1 "" 1.2 

2,1 ^ 2,2 ^ 2,3 

^ 3,2 '^ 3,3 '^ 3,4 


n-l,n-2 n-l,n-l n-l,n 

z z 

n,n-l n,n 



bA Lb J 

K.s] 1 K.e] 
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The sign convention used in this paper to obtain the mass 
matrix [Z] is shown in Figure 2* 



The beam mass matrix obtained in this paper is applicable 
for either the pitch or yaw plane with a change of variables and 
possible sign changes. For the axis system shown in Figure 3, 
the following variables would be used. Note that in order to 
have a right-hand system such as that shown in Figure 3, the sign 
of the and |^Z^^ partitions of [Z] from this subroutine 

would have to be changed for the yaw plane. 


This 
Paper 
(A, X, 9) 

Axis System 

of Figure 3 
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(z. X, Gy) 

Yaw 

(>'• *• “z) 

X 

X 

X 

6 

6 

6 


z 

y 

0 

e 
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y 

z 

m(x) 

m(x) 
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r(x) 

r (x) 
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y 

z 

M 

M 


c 

c 

c 
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iCg 

iCg 

c 

c,y- 
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e 

Z 



The beam’s total mass, center of gravity, and moment of iner- 
tia about the beam’s center of gravity are also calculated and 
printed. 


DESCRIPTION OF TECHNIQUE 


The replacement of distributed mass, distributed rotary iner- 
tia, and concentrated mass items of a beam by a mass matrix is 
obtained using a kinetic energy approach as follows. 


For small deflections normal to the longitudinal body axis (x) 
of the beam shown in Figure 2, the kinetic energy is defined by 


T = 


h I l^ni(x) 6^(x,t) + r(x) 6^(x,t)j 


dx + 


V fw e2(t) 

/ ^ [ c c c c 


] 


( 1 ) 


where 


m(x) is the distributed mass, 

r(x) is th 'istributed rotary inertia, 

M is the mass of a concentrated item, 
c 

C2 

I ^ is the rotary inertia of a concentrated item about 
its own center of gravity, 

* 

6 is the time rate of change of elastic deflection 
normal to the bocy x-axls, and is referred to as 
lateral velocity in this paper, 
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0 is the angular velocity about an axis perpendicular to 
the paper, 

A is the velocity of a concentrated item and will be ex- 
plained later, 

X is the undeformed longitudinal axis of the beam, 
t is time, 

Xg is the starting x-station of the beam, and 

X- is the ending x-station of the beam. 
t» 

The finite summation is over the number of concentrated mass items. 

The technique used here to integrate Equation (1) is essen- 
tially a modification of the Rayleigh-Ritz method. In the basic 
Rayleigh-Ritz method the elastic deflection of a beam is repre- 
sented by a summation of assumed displacement functions each 
weighted by a generalized coordinate representing the contribu- 
tion of each of the assumed functions. Each displacement func- 
tion is assumed over the entire beam length and considerable skill 
is required to choose these functions. The modification used here 
is to assume a simple displacement function between consecutive 
panel points and to use the panel point displacements as the gen- 
eralized coordinates. A cubic displacement function will be as- 
sumed here. This is illustrated in Figure 

are consecutive panel point stations. The 
points k and k+1 is referred to as bay k. 


4 where x, and x, . - 
k k+1 

region between panel 


6(x) 





X 


k+1 


X 


Figure 4 Cubic Displacement Function 
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The cubic displacement function in bay k is determined as fol- 
lows. The lateral displacement is given by 

or 


6(x) = A|x - + b|x - Xj^j^ + c| 


6 (H) = aH ^ + bH ^ 4* cH + d 


where 


with 


= r2 H 



H = 



\ “ \+l ■ 


( 2 ) 

(2a) 


C2b) 


U is a local nor.dlmensional bay coordinate to ale in the solution 
of Equation (1) . The angular displacement Is obtained as the 
geometric derivative of the lateral displacement. That is* 


or 


G(x) 


96(x) 

3x 


0(H) 


d 6(H) 

■ L. dH 
k 


= :^[-3h2 


a 


b 


-2H -1 01 


c 

d 


(3) 
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The coefficients a, b, c, d are determined as follows from Equa- 
tions ( 2 ) and (3) . 



L 1 0 0 0 J 


Thus the lateral and angular displacement in bay k is determined 
in terms of the adjacent panel point lateral and angular displace-- 
ments* Similarly, the lateral and angular velocity in bay k is 
given by 


6(H,t) = [r3 r2 H 


1 ] W 



(A) 
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and 


0(H, t) = ^ [-3 h2 


-2H -1 0] [«i-] 


6,(t) 

6k(0 


(5) 


The previous expression for the angular velocity is slightly in 
error because a rotation due to shear deformation is included. 

A knowledge of the stiffness distribution would be necessary to 
determine the relative contribution of bending rotation and shear 
deformation to the angular velocity. This error will be small 
for beam- type structures and will be ignored. 


The kinetic energy of distributed mass, distributed rotary 
inertia, and concentrated mass items will be considered separately. 
The distributed mass, m(x) , is considered first. The geometry for 
a line segment of distributed mass is shown in Figure 5. 


m(x) 



*1 X2 

Figure 5 Line Segment Geometry 


The equation for a straight line segment as 

5 Is 

m(x) = mi + (x - xi) (m 2 - mi)/(x 2 - 


shown In Figure 
xi) , 


or 


m(H) = WH + mi - WHi 


( 6 ) 


by using the nondlmenslonal coordinate of Equation (2a) and 

W = (m 2 " mi)/(H2 - Hi). 


(7) 
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Modifying Equation (1) to use the nondlmensional coordinate of 
Equation (2a) and substituting Equations (4) and (6) gives the 
kinetic energy of the distributed mass represented by one line 
segment 1 In bay k as 


J,(o 

/ 

H 




M* 

\ 



ivi^L 

f' 

h5 




WH + m - Wll <1 h1i. ) 


0,(0 

^ 1 k 




U'"' 

H 

1 ” ”1 

■•.,(o 

0,,.j^(t) 

V 

P 

h3 


H 

1 

/ 



The subscripts p and q have been Introduced to handle the possl*~ 
blllty of a line segment extending past the bay limits. Thus, 

X (or H \ is the greater of x, or , and x /or H \ is the lesser 

p I p/ Ik q \ q/ 

of x^ or X, . , . Similarly, m is either m, or m, and m is either 
2 k+1 P Ik q 

m2 or integration is continued for the line segment in 

adjacent bays, if necessary until the entire line segment has 
been used. Performing the integration of Equation (8) yields 


* ■‘I* 


6,(0 




9^(0 


~ 


" 

, k k+1 ^k , k+n ^k , k+n+1 



^k+l.k+1 \+l.k+n ^k+1, k+n+1 



^k+n,k+n ^k+n, k+n+1 


0,(0 

^k+n+1, k+n+1 


k+1 



(9) 


where 


\,k 
^k , k+1 
^k,k+n 
\,k+n+l 
\+l,k+l 
\+l,k+n 
\+l,k+n+l 


Fi - 6F3 + 4F4 + Pj 

3F3 - 2F4 - Pi 

(-F2 + 2F3 - F4 - P2) \ 

(F3 - F4 - P3) 

Pi 

"3 H 


(9a) 

(9b) 

(9c) 

(9d) 

(9e) 

(9f) 

(9g) 
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*lt+n,k+n 

*k+n,k+n+l 

\+n+l,k+r+l 

Pi 


(2Fb - 4F6 + F 7 + Fa - AFg + AFiq) 
(Fa - 3Fe + Fy - Fg + 2 Fiq) 

(-2Fe + F 7 + Fio) \ 

-12Fe + 4Fy + 9Fio 


?2 * 2 F 5 - 7Fe + 2 F 7 - 3Fg + 6 F 10 


P 3 » -5Fe + 2 F 7 + 3Fio 


fa - E 3 

Fg = Lj^ E 4 
PlO ® ^ ^5 

Ej - «(»r^ - * 

("p - ““p) H - “J)/^ 

%) (\ - «p) 

“p " ,(‘p ' “eI/'x 
■ (X, - \)/hc 

'•k ■ Xk+l • *X- X““ 
n Is the number of panel points • 


(9h) 

(9i) 

(9j) 

(9k) 

(9t) 

C9m) 

(9n) 

(9p) 

(9q) 

(9r) 

(9s) 

(9t) 

(9u) 

(9v) 

(9w) 


It is possible to simplify these equations for z since Fg = F 3 , 

Fg ■ F 4 , and Fjo • F 5 . This was not done so that these same equa- 
tions for z could also be used for distributed rotary inertia and 
concentrated mass items and thus save coding instructions in the 
subroutine. The kernel matrix in the triple matrix product of 
Equation (9) is the mass matrix which replaces the distributed 
mass represented by one line segment i in bay k. 
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The kinetic energy of Che distributed rotary inertia, r(x), 
is considered next. The geometry for one rotary inertia line seg- 
ment is similar to Chat given in Figure 5 for distributed mass. 
Similarly, 

r(H) » WH + ri - WHi (10) 

using the nondimensional coordinate of Equation (2a) and 

W - (r2 - ri)/(H2 - Hi). (11) 

Modifying Equation (1) to use the nondimensional coordinate of 
Equation (2a) and substituting Equations (5) and (10) gives the 
kinetic energy of the distributed rotary inertia represented by 
one line segment i in bav k as 



The subscripts p and q have similar meanings as has been previ- 
ously discussed* for distributed mass. Performing the integration 
of Equation (12) yields results identical to Equations (9) thru 
(9m) and similar to Equations (9s) thru (9w) as obtained previ- 
ously for distributed mass* These equations will not be repeated 
here* The differences from the distributed mass equations are: 


= 0 (j = 1, 4) 

(13a) 

Fs = 3E3/Lj^ 

(13b) 

Fe = eE^jL^ 

(13c) 

F7 = 

(13d) 

Fe = El/Ll, 

(13e) 

Fg = 2 E 2 j'L^i and 

(I3f) 

FlO “ 

(13g) 


The kernel matrix in the triple matrix product of Equation (9) is 
then the mass matrix which replaces the distributed rotary inertia 
represented by one line segment 1 in bay k. 
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The concentrated mass items are considered last. These items 
require some investigation because the item may be connected to 
the beam at an attach point a by an arm (assumed rigid) as shown 
in Figure 6. 


6 



The square of the velocity of concentrated item c is 

A2(t) = |6^(t) - 0^(t)j2 + |ne^(t)j2 

From Equation (1) , using the velocity expression of Equation 
(14 ) 9 the kinetic energy for a concentrated item c becomes: 


(14) 


'c “ ''K - "a) “a<"> ^[(*0 - *a) 


2 + d2 


02 (t) + 1*^® e^(t)> 
a c a 7 

(15) 


where use has been made o£ Che relationship Chat 0 (c) = 0 (c) . 

c a 

In further work, the offset D will be assumed to be zero. The 
reason for this is that any mass offset from the beam x-axis means 
that the longitudinal principal axis of the beam does not coincide 
with the beam x-axis. A coupled transverse-longitudinal analysis 
wou?.d be required to accurately describe the kinetic energy of the 
beam for D 0. However, if the analyst wishes to include the ef- 
fect of an offset D, inspection of Equation C15) shows that 

+ M could be used. With the above conditions, the kinetic 
c c 

energy for a concentrated item can be expressed in matrix notation 
as 


T 

c 


a 


0a(t) 



-M 



-H lx. - X \ 
c\ c aj 


M/x - x\ ^ + I 


/X “ X \ 

c^ c aj 


eg 



(t)i 


a 


0 (t) 


a 


(16) 
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The lateral and rotational velocities at the attach point a are 
given in terms of the velocities at adjacent panel points k and 
k+1 from Equations (4) and (5), respectively* Using these equa- 
tions and modifying Equation (16) to use the nondimensional co- 
ordinate of Equation (2a) , the kinetic energy of one concentrated 
item c in bay k is identical tw Equations (9) thru (9m) as ob- 
tained previously for distributed mass* These equations will not 
be repeated here* The differences from the distributed mass 
equations are: 


Fi = M 
^ c 


(17a) 

F2 = M H 
^ c c 


(17b) 

Fa = H M 
^ a c 

1 



(17c) 

Fi* = M 

(3H - 2H \ 

(17d) 

^ a c 

V c a/ 



Fs = h2 {m H /3H - 2H \ + 

^ a[cc\c a/ c/ic 

= “a|«c(“c - “4" * "c'/'il 
Fs = 

F, = H^|m^ - H,) + 2l‘yL^| 

F,„ - H||m^( 2H^ - * «=S/L2j 



C17h) 
(171) 
(17j) 
. (17k) 

(171) 


(17e) 

(17f) 

(17g) 


C17m) 
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The mass matrix for the entire beam Is obtained by evaluating 
£*i.atlons C9) thru (9ir) for each line segment of distributed mass* 
£q lations (9) thru (9m) > (9s) thru (9v) , and (13a) thru (13g) for 
each line segment of distributed rotary inertia » and Equations 
(9, thru (9m) and Cl7a) thru (17m) for each concentrated mass 
ixcm« Every resulting 4x4 bay mass matrix is added to previous 
4x4 bay mass matrices at like panel points to form the mass matrix 
for the entire beam. 


The total mass properties of the beam are calculated by the 
fol Lowing triple matrix product. 



vhere 


( 1 } 

{ 0 } 

{PP} 



is a column of ones* 
is a column of zeroes, 

is a column of the panel point x-stations, 
are partitions of the mass matrix [Z], 

is the beam mass. 


P]^ is the first moment of the beam about x - 0, and ' 
is the moment of inertia of the beam about x = 0. 


From this data, the center of gravity of the beam is calculated 
from 


X 

eg 



ana the moment of inertia of the beam about x is calculated 



X 


2 

eg* 


om 
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SPECIAL CA SES 

Two cases are worthy of mention. In the first case, a 
lirutco'hj varying segment of distributed mass extending to the 


bay limits | 

[*k 

of: 


“0\ * ’^1 



72\ + 240n^+l \ ~ + ^^'\+l)^k ^°\+l) 

1 (5*\ ^”k+l)^k -(H H+l)^ 

(sym) I 

1 P’V'" 5\+l)^k 


where 


In the second case, a uniform |i.e., = *\+l "" °^) segment of 

distributed mass extending to the bay limits and has- 

a mass matrix of ’ ^ ^ 

TlSe 54 1 -22L, 13L, 1 


156 I -13Lj^ 22Lj^ 


420 1 

I K 

(sym) I 


4l2 -3l2 

k k 


where 


* ^1 -LI ” ^ • 

k k+1 K 

This last result agrees with Archer’s consistent mass matrix for 
a uniform (mass and stiffness) beam segment given in Reference 1* 
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Subroutine MASS2A takes distributed mass of fluid in a con- 
tainer and the slosh mass of the same fluid and replaces them with 
a mass matrix. The mass matrix includes coupling between these 
two types of masses. The elements in the mass matrix for the dis- 
tributed mass and the coupl^n 3 terms are representative inertial 
value's at selected points on the beam. These elements are cal- 
culated by assuming (1) a cubic function between pairs of the 
selected points to describe the container’s lateral elastic de- 
flection> and (2) a uniform motion of the fluid slosh relative to 
the container’s deflected longitudinal axis. 

The X stations of the selected points (panel points) are given 
in {PP}. These x stations must be in increasing order. 

The distributed mass» m(x) » is assumed to be piecewise linear 
and is represented by straight line segments as shown in Figure 1. 



The X stations of the end points for the line segments giving the 
distributed mass are independent of the panel point x stations. 
However, the distributed mass must be within the panel point 
limits. The line segments representing the distributed mass may 
or may not be Joined; however, there must not be any x voids or 
overlaps. The distributed mass is defined in [DMASS]. Each row 
of [DMASS] represents one nonvertical line segment. The form of 
each row of [DMASS] is [xj X 2 ^ 2 ] where xj, give the first 
end point and x^* m 2 give the second end point of a line segment. 
The X 2 of row 1 of [DMASS] must be equal to of row 2, etc. 

To aid the analyst, the fluid level (FLEVEL) is input so that 
[DMASS] can describe the complete container fluid. Only the 
fluid from the fluid level to the container bottom will be used 
in calculating the mass matrix. Another feature is that [DMASS] 
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can be used to define the cross-sectional area of the complete 
container and the input scalar CONVRT used to define the fluid 
density. 

It is assumed that there is no distrihutei rotary Inertia for 
a fluid as there is for a structure. 

The slosh mass (SMASS) of the fluid is a value that is sepa- 
rately obtained either from test or mathematical modeling tech- 
niques . 


The calculated inertial terms are placed in a mass matrix 
given by [Z]. The form of [Z] is 


^^^2n+l,2n+l 





h.e] 



pa..] 





. i'l-i 


"i 

ps.e) 

T ! 

i : 

Z 

! 


where n is the number of panel points. Matrix [Z] is symmetric, 

=s The partition form of [Z] results from the two 

generalized coordinates (lateral translation 6 and rotation 8 ) at 
each panel point arranged with all translation coordinates first 
followed by all rotation coordinates. These panel point coor- 
deinates are then followed by the single fluid slosh coordinate. 
Each panel point partition of [Z] is square and tri-diagonal, for 
example. 
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The element Z Is the input slosh mass (SMASS) . 

S ) s 

The sign convention used ir this paper to obtain the mass 
matrix [Z] is shown in Figure 2. 



Figure 2 Sign Convention for This Paper 

The mass matrix obtained in this paper is applicable for either 
the pitch or yaw plane with a change of variables and possible 
sign changes. For the axis system shown in Figure 3, the follow- 
ing variables would be used. 


This 
Paper 
(6 ,x,0) 

Axis System of 
Figure 3 

Pitch 

Yaw 

(y.x.ej 

X 

X 

X 

6 

6 

6 


z 

y 

0 

6 

-e 


y 

z 

m(x) 

m(x) 

m(x) 
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t 


0 

z 



^0 

X 


Figure 3 Beam Axis System (Right Hand) 

Note that in order to have a right-hand system such as that shown 
in Figure 3, the sign of the and |z^^gjT 

partitions of [Z] from this subroutine would have to be changed 
for the yaw plane. 

The fluid's total mass, center of gravity, and moment of 
inertia about the fluid's center of gravity are also calculated 
and printed. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed mass and the slosh mass of 
fluid in a container by a mass matrix is obtained using a kinetic 
energy approach as follows. 

For small deflections normal to the longitudinal body axis (x) 
of a container, the kinetic energy of the fluid in the container 
shown in Figure 2 is defined by 


T 



6^(x,t) + 


( 


6 (x,t) + 6 (x,t) 
c s 



dx 


( 1 ) 
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where 


m^(x) is the distributed inert mass of the fluid, 

m (x) is the distributed slosh mass of the fluid, 
s * 

6 (x,t) is the time rate of change of the elastic deflec- 
tion of the container's longitudinal axis normal 
to the body x-axis, 

• 

6 (x,t) is the time rate of change of the fluid sloshing 
motion relative to the container's deflected 
longitudinal axis, 

X is the undeformed longitudinal axis of the con- 
tainer , 

t is time, 

X is the starting x-station of the fluid, and 

4 > 

is the ending x-statlon of the fluid. 


The distributed mass, m(x). Is the sum of the distributed Inert 
mass and the distributed slosh mass, that Is, 


Define 


ra(x) =» m (x) + m (x) . 
1 s 

SMASS 

^OT 


where SMASS 

and Is 

rO( tine* 


is the slosh mass and is input to the subroutine, 
the fluid total mass and is calculated in the sub- 


Assume that the distributed slosh mass and distributed mass vary 
in the same ratio as the total respective masses, that is, 


m(x) 


= R 
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Then Equation (1) nan be written as 


T 



- R) 6^(x,t) + R 


^<S^(x,t) 






dx. 


(la) 


The technique used here to describe the elastic defl^c^ion 
the container is essentially a modification of the Rayleigh-Ri< 
method. In the basic Rayleigh-Ritz method the elastic deflect ic 
of the container is represented by a summation of assumed displace- 
ment functions each weighted by a generalized coordinate represent- 
ing the contribution of each of the assumed functions . Each dis- 
placement function is assumed over the entire container length 
and considerable skill is required to choose these functions. The 
modification used here is to assume a simple displacement func- 
tion between consecutive panel points and to use the panel point 
displacements as the generalized coordinates. A ovihio displace- 
ment function will be assumed here. This is illustrated in Figure 
4 where and consecutive panel point stations. The 

region between panel points k and k+1 is referred to as bay k« 



Figure 4 Cubic Displacement Function 

The cubic displacement function in bay k is determined as follows. 
The lateral displacement is given by 
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<S (x) 
c 


= A(x - - Xj^)2 + C(x - + 


or 




“aH^+bH^ + ci 


[h3 h2 H 1] 


a 

b 

c 

d 


where 


“ “I* ■ \)/‘ic 


with 


L, =* X, 


k+1 


- X, 


k* 


H is a local nondimenslonal bay coordinate to aid In the solu~ 
tlon of Equation (1). The angular displacement is obtained as 
the geometric derivative of the lateral displacement. That is, 

'd 6 (x) 

d— 


or 




d 6 (U) 
c 

L, dK 
k 


^ [-3h 2 -2U -i 0] 


1 
I b 

c 


lAi 


( 2 ) 

(2a) 


(2b) 


(3) 
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The coefficients a, b, c, d are determined froni Equations (2) and 
(3) as follows: 



Thus, the lateral displacement in bay k is determined in terms of 
the adjacent panel point lateral and angular displacements* Sim- 
ilarly, the lateral velocity in bay k is given by 




[h3 H 1] [4-] 



(4) 
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The fluid sloshing motion, ^^(x,t), relative to the container's 

deflected longitudinal axis is given as the product of a displace- 
ment function, f^(x), and the slosh amplitude, A^(t). A unifovm 

displacement function, f (x) » G, will be assumed here. The value 

s 

of G will be determined later. The time derivative of the slosh- 
ing motion is then 

• • 

6 (x,t) = G A^(t). (5) 

s s 


Combining Equations (4) and (5) gives the total velocity of 
the fluid as 


6g(H,t) 


+ 6g(H,t) = 


[H^ H 1 G] [4;] 


6^(t) 






where 


M * 


2 

-3 


0 


1 

0 


-‘•k 0 

^ ^‘■k '•k 0 

0 -L. 0 0 

k 

0 0 0 0 

0 0 0 1 


( 6 ) 


The distributed mass of the fluid will be considered next. 

The geometry for a line segment of distributed mass is shown below 
in Figure 5. 
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Figure 5 Line Segment Geometry 

The equation for a straight line segment as shown in Figure 5 is 
m(x) * mx + (x - xi) (m 2 - mx)/(x 2 - xi) 


or 


m(H) * WH + mi - WHi (7) 

by using the nondimensional coordinate of Equation (2a) and 

W - (m 2 - mi)/(H2 - Hi). (8) 

Modifying Equation (la) to use the nondimens ional coordinate of 
Equation (2a) and substituting Equations (4), (6), and (7) ^ives 
the kinetic energy of the distributed mass of fluid represented 
by one line segment i in bay k as 


i^(t) 





H'* 


H^GR 

|WH + nip - WHpJdH^ 

[c] 




Jh 

p 


11“ 



H'’gK 




o^(t) 



u'* 

h3 

H2 

U 

HGR 



Vt) 



i 


H2 

U 

1 

GR 





\ 

^H^GR 

H^GR 

HGR 

GR 

g2r _ 

1 


A (t) 

L 0 J 


( 9 ) 
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The subscripts p and q have been introduced to handle the possi- 
bility of a line segment extending past the bay limits. Thus, 

Xp ^or is the greater of x^^ or Xj^, and x^ |or is the lesser 

of X 2 or Similarly, m^ is either m^^ or nij^, and m^ is either 

m^ or The integration is continued for the line segment in 

adjacent bays, if necessary, until the entire line segment has been 
used. Performing the integration of Equation (9) yields 


where 


e^<t) 


! 7 , 


*k,k+i *k,k^^l ^,2n+l 

*k+l*k+l ^+I,2n+1 

\*n,k*n. *i^,k-Ht4-i *k4n,2xH-l 

^♦trU.kfnVl %:+Tt4^l,2a^l| 






k»k 

^,k+l 

^k,k+n 

\.k+n+l 

*k+l,k+l 

*k+l,k+n 

\+l,k+iri-l 

*k+n,k+n 

\+n,k+n+l 

*k+n+l,k+n+l 

\, 2 n+l 

*k+l, 2 n+l 


= - 6 F 3 + 4F^ + P 3 

= 3 F 3 - 2 F^ - P 3 

* ^'^3 - '4 - ^2)\ 
■ ^ - ^3)Hc 


("3 


= P 


1 


("3- 

(-"4 


4 3K5 - 3Fj 4 F,)l2 


- ^”6*‘'K 

- 3 F 3 + 2F^|GR 

- 2 F 4 )GR 


(10) 


( 10 a) 

( 10 b) 

( 10 c) 

(lOd) 

(lOe) 

(lOf) 

(lOg) 

(lOh) 

(lOi) 

(lOj) 

(lOk) 

(lOi) 
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*k+n,2n+l 

- (-'^2 + ^3 - ^4)^ ® 

(10m) 

*k+n+l,2n+l 

■ (^3 - “ 

(lOn) 

*2n+l,2n+l 

= G^R 

(lOp) 


= 9Fj -12Fg + 4Fy 

(lOq) 

"2 

= -3F, + 8F, - 7F, + 2F- 
4 5 o 7 

(lOr) 

"3 

- - 5^6 ^■'7 

(10s) 

F. 

J 

' V 



/m - WH \ /h^ - H^)/j|L. 

(lOt) 


\ p p) \ q p// / ^ 


W 

\ q p) /I q p) 

(lOu) 

u 

p 

■ (-P - »k)/‘-k 

(lOv) 

H 

q 

' (% -*k)/\ 

(lOw) 


' \n - *k- 

(lOx) 


n is the number of panel points* 


G is determined as follows. Because ^ _ is defined to be 

2n+l,^n+l 

equal to the slosh mass» SMASS, 



where the summation is over the number of segments giving the dis- 
tributed mass of the fluid. From Equation (lOt), 

■ i (», - *p) C'p ^ %} 

which is the fluid mass in the inteval x to x . The summation 

p q 

of is then the fluid total mass, ^ previously de- 

fined as the ratio SMASS/M^^^. Therefore 


G * 1. 


( 11 ) 
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The kernel matrix in the triple matrix product of Equation (10) 
Is the mass matrix that replaces the distributed mass of fluid 
represented by one line segment 1 in bay k. 

The mass matrix for the entire fluid is obtained by evaluating 
Equations (10) for each line segment of distributed mass. Every 
resulting 3x5 bay mass matrix is added to previous 5x5 bay mass 
matrices at like panel points to form the mass matrix for the 
entire fluid. 

The technique just described was developed by Mr. Carl Bodley 
and Mr. Herbert Wilkening. 


The total mass properties of the fluid are calculated by the 
following triple matrix product: 


U) 

{PP} 

I 

>..] j 



{1} 1 {PP} 

1 


«T 

1 

{0} ‘ 

1 

{-1} 


J^Oi] 1 

w. 


1 

(0) j {-1} 




where 


{1} is a column of ones, 

{0} is a column of zeroes, 

{FP} is a column of the panel point x stations, 
etc., are partition of the mass matrix [Z], 

M,^, is the fluid mass, 
o 

is the first moment of the luid about x = 0, and 
o 

is the moment of inertia of the fluid about x = 0. 


From this data 


and the moment 


the center of gravity of the fluid 


eg 


■ ‘•?A 


of inertia of the fluid about x 


eg 


= l" - M., x^ . 

i 1 T eg 


is calculated from 


is caicul.iLed from 
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Subroutine MODEl calculates the mode shapes and natural fre** 
quencles of a structure using its mass and stiffness matrices. 

The method of Jacobi is used in the solution of the general 
eigenvalue problem 

^ ^NxN 

where the mass matrix is given by [A], the stiffness matrix by 
ISj, and the size of the matrices by N. The mass and stiffness 
matrices must both be real and symmetric. The mass matrix must 
also be positive definite. The i,i element of the calculated 

th 

diagonal matrix T J is the i circular frequency squared. 

th 

The corresponding mode shape is the i column of ['t']. Because 
Jacobi's method is used, all mode shapes and frequencies are cal- 
culated. The circular frequency squared, , the circular fre- 
quency, i*j, and the frequency, f = uj/Ztt, are output from the sub- 
routine in increasing order of magnitude in {W2} , {W}, and (FREQ;, 
respectively. The mode shapes are altered so that the first ele- 
ment in each column is positive and normalized (automatically due to 

T 

Jacobi's method) so that [A] Orthogonality checks 

T T 

[‘‘r'j [A] and ['I’] [S] [P] are calculated and a summary of 

these results is printed (on option) • 


DESCRIPTION OF TECHNIQUE 

Consider a discrete coordinate model of a structural system 
having N degrees of freedom. The undamped, free equations of 
motion for small vibrations about a position of stable equilibrium 
can be expressed in matrix notation as 


(h(t)L 


{h(t)}„„, = {Oh 


As before, [A] is the mass matrix and [S] is the stiffness matrix. 
(h(t)} is a column of discrete coordinate displacements. The 
time dependence in Equation (2) can be removed by the transforma- 
tion {h(t)i = e^^^ because the solutions of Equation (2) are 

assumed to be harmonic in time. Thus, for a nontrivial solution, 
Equation (2) becomes 


(-'.Ma] + [S]) (v) - (0}. 


(3) 
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[A] { 


Equation (3) is recognized as a general eigenvalue problem of 
order N with eigenvector (mode shape) and eigenvalue w* 
(circular frequency squared). There are N values of u' and {.} 
so that the expansion of Equation (3) to include all solutions 
can be expressed as 


I?}! 


I I 

(4)}2 I I 


2 

(U2 



or 


[A] [$] = 

Premultiplying by [A]“^ gives 

IA]-1 tS] 14-] = [<t>] 


= IS] 

{0>1 




(4) 

in- 


(Aa) 

ro.2 j 


(5) 


as the eigenvalue problem to be solved in this paper. Jacobi's 
method cannot be used directly on Equation (5) because even 
though [A] and [S] are symmetric, their product in general is 
not. 


The first operation in the solution of Equation (5) is to 
decompose (reference subroutine DCOMl) the mass matrix into tri- 
angular factors, that is. 


[A] = [U]" [U] 


where [U] is an upper ^ triangular matrix. Using [A]"^ * ([U]^ 

[U])“^ = [U]“^ premultiplying Equation [5] by [U] 

gives ' ^ 


( 6 ) 



[s] in 


[u] in 


(7) 
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Define 


[t] = [U] [4>] (8) 

from which 

W = [U]-l [?]. (9) 

Using the fact that ([U]^) = we get 

[D] [T] = [T] (10) 

where 

[Dl - ([U]-l)^ [S] [Ul-1 (10a) 

and la rafarrad to aa tha dynamic matrix In the aubroutlna. Be- 
cauaa the atiffnaaa matrix [S] la aymmetrlc, [D] la aymmatrlc. 

Equation (10) is the eigenvalue problem from which the eigenvalues 
r (A)^ J and eigenvectors [$] are calculated by Jacobies method 
(reference Subroutine EIGNl). The product of equation (9) to 
obtain the mode shapes [i] is accomplished by using [uf as the 
Initial eigenvectors In Subroutine EIGNl ^ thus eliminating a 
matrix multiplication. 

Two orthogonality checks of the mode shapes are calculated in 
the subroutine aa follows. Firsts a sunanary of the results of 

T 

[4»1 [A] are printed. This calculation checks the orthogonal- 

ity of the mode shapes on the mass matrix 

and the normalization |i.e., [A] 

of [4»]^ [S] [4»] are compared with J . (They should be equal.) 

T 

This check results from premultiplying Equation (4a) by ['1'] 

T 

and assuming [^] [A] [^] = [I]. Only the upper half of [A] and 

[S] are used in the calculation of the mode shapes and frequencies 
but all of [A] and [S] are used in the orthogonality checks. Thus, 
errors in the checks can result from nonsymmetric [A] or [S]. 


|i.e. , (4'}^ [A] {y;^ = oj 
= l). Next, the results 
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Subroutine MODEIA calculates the mode shapes and natural fre- 
quencies of a structure using its mass and stiffness matrices. 

The method of Jacobi is used in the solution of the general 
eigenvalue problem 




NxN 


diagonal matrix 
squared 


1 deiinite. 

r= + “d 


where the mass matrix is given by [A], the stiffness matrix by 
(S], the size of the matrices by N, and c is a constant automat- 
ically calculated in the subroutine. The mass and stiffness 
matrices must both be real and symmetric. The mass matrix must 
also be positive definite. The i,i element of the calculated 

contains the i^^ circular frequency 

th 

The corresponding mode shape is the i column of [v]. 
Because Jacobi’s method is used, all mode shapes and frequencies 
are calculated. The circular frequency squared, , the circular 
frequency, u, and the frequency, f = a)/2TT, are output from the 
subroutine in increasing order of magnitude in {W2}, {W}, and 
{FREQ}, respectively. The mode shapes are altered so that the 
first element in each column is positive and normalized (auto- 

T 

matically due to Jacobi’s method) so that [A] = 

T T 

Orthcgonali ty checks [4>] {A] [1^] and ['J'] 

and a summary of these results is printed 


i 

[S] [i^] are calculated 
(on option) « 


( 1 ) 


DESCRIPTION OF TECHNIQUE 


Consider a discrete coordinate model of a structural system 
having N degrees of freedom. The undamped, free equations of 
motion for small vibrations about a position of stable equilibrium 
can be expressed in matrix notation as 






= {0} 


Nxl‘ 


( 2 ) 


As before, [A] is the mass matrix and [S] is the stiffness matrix. 
(h(t)} is a column of discrete coordinate displacements. The time 
dependence in Equation (2) can be removed by the transformation 

{h(t)} = a{(J)} e"^^^ because the solutions of Equation (2) are as- 
sumed to be harmonic in time. Thus, for a nontrivial solution, 
Equation (2) becomes 


(-oj^[A] + [S]) {.*,} = {0}. 


(3) 
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Equation (3) is recognized as a general eigenvalue problem of 
order N with eigenvector (mode shape) and eigenvalue (cir- 

cular frequency squared). There are ll values of -jj" and so 
that t.ie expansion of Equation (3) to include all solutions can 
be expressed as 


[A] 


i 1 ! ' 


I V V ; , 


OJ"’ 


- [SI 


(4^}] 1 \ “Mi. 



(4) 


or 


[A] [c] r j = [s] [1^]. 


(Aa) 


Adding c[A] [y], where c is a constant, to both sides of the 
equation results in 


(c[A] + [S])-l [A] = ['H 



c + wg 


(5) 


as the eigenvalue problem to be solved in this paper. The con- 
stant, c, is calculated from 


_ Norm (S) 

Norm (A) 

where the Norm of a matrix is defined as the sum of the absolute 
values of the elements of the matrix. Jacobi's method cannot be 
used directly on Equation (5) because even though [A] and [S] are 
symmetric, their product in ge jral is not. 

The first operation in the solution of Equation (5) is to de- 
compose (reference Subroutine UCOMl) into triangular faLtors, 
that is, 


C [A] + IS] - [U]^ [U] 
where [U] is an upper triangular matrix. 

Uditig (c (A] + [S]) ^ [U] ^ ([U]^) premultiplying Equation (5) 

T -1 

and using a superscript on t<j)] gives v[U]‘) [A] l'|>] “ 


A 


r 1 

C+0) 


( 7 ) 
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Next, we define 

m =■ [U] [4] 

from which 

[^] = [u]-i [T]. 

Using the fact au]”^)"^ - we get 

[D] [j] = [T] r A J 

where 

[.-] - ([u]"^)”^ [A] [U]"^ 

and 

M = 


C + 


( 8 ) 

(9) 


( 10 ) 


(10a) 


(10b) 


[D] la r«f«rr«d to as the dynamic matrix In the eubroutln«s Be- 
cause the mass matrix [A] Is symmetric, [Dl Is symmetriCe Equa- 
tion (10) Is the eigenvalue problem from which the eigenvalues 
and eigenvectors [<t] are calculated by Jacobi's method 
(reference Subroutine EIGNl) . The product of Equation (9) to 
obtain the mode shapes [$] Is accomplished by usln^ [bl"l as the 
Initial eigenvectors In Subroutine EIGNl, thus eliminating a 
matrix multiplication, Tiie mode shapes [^] must be renormalized 

T X 

because where [<t>] [A] [<})] ■ [I] is wanted, we have [<{)] [A] [<{)]■ 

[X], This last telation is obtained from Equations (10), (10a), 

and (S) and using the orthon.')rmal properties of [(J)],, that is, 

^ -1 ^ T 

■ 14^] • The final renormalized mode i is obtained from 


(11) 


( 12 ) 

Two orthogonality checks of the mode shapes are calculated 
(on option) in the subroutine as follows, First a su>:iiary of the 
T 

results of (A] [(t>] are printed. This calculation checks 


«), 

• 

th 2 

The circulai' frequency squared, to is obtained from Equa- 
tion (10b) as ^ ^ 

^ - 1;; “ 
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the orthogonality of the mode shapes on the mass matrix (i.e. 

[A] {<)}j - 0) and the normalization (i.e. [A) «1) . 

Nextf the results of 14>j'‘ [S] [4>] are compared with 

(They should be eqifal.) This check results from premultiplying 

T T 

Equation (4a) by [^] and assuming [^] [A] [<1>] « [1]. Only the 
upper half of [A] and [Sj are used in the calculation of the node 
shapes and frequencies* but all of [A] and [S] are used in the 
orthogonality checks. Thus* errors in the checks can result from 
nonsymmetric [Al or [S]. 
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Subroutine MODEIB calculates the mode shapes and natural fre- 
quencies of a structure using its mass and flexibility matrices. 
The method of Jacobi is used in the solution of the general eiger-' 
value problem 


tAl 




[^1 





( 1 ) 


where the mass matrix is given by [A], the flexibility matrix by 
[El, and the size of the matrices by N. The mass and flexibility 
matrices must both be real and symmetric. The mass matrix must 
also be positive definite. The i,i element of the calculated 

diagonal matrix is the reciprocal of the i^^ circular fre- 

Lw J th 

quency squared. The corresponding mode shape is the i column 

of [^]* Because Jacobi's method is used, all mode shapes and 
frequencies are calcuxated. The circular frequency squared, 
the circular frequency, co, and the frequency, f - oj/2r, are output 
from the subroutine in increasing order of magnitude in \W2}, 

(V), and {FREQ}, respectively. Rigid body frequencies and mode 
shapes «/ill be in the lact positions of {W2}, (W), {FREQ}, and 

[^]. This is because ~ will be calculated as zero (or a small 

number due to computer roundoff error) corresponding to the rigid 
body modes. The mode shapes are altered so chat the first ele- 
ment in each column is positive and normalized (automatically due 

T 

to Jacobi 3 method) so that (4»}^ [A] {v)^ - 1. Orthogonality 

checks [AJ and (4']^ {A] [E] [A] [4>] are calculated and 

a summary of these results is printed (on option). 


DESCRIPTION OF T£CHi:iQU£ 


Consider a discrete coordinate model of a structural system 
having N degrees of freedom. The undamped, free equations of mo- 
tion for Fmall vibrations about a position of stable equilil'^riu.n 
can be expressed in matrix notation as 


As before, [A] IS the mass matrix, [S] is the stiffness matrix, 
and {h(t)} is a cclumn of discrete coordinate displacements « The 
time dependence in Equation (2) can removed by the transforma- 
tion {h(t)) * oil';)} because the solutions of Equation (2) are 

assumed to be harmonic in time. Thus, for a nontrivial solution, 
Equation (2) becomes 
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(-to^lAj + [Si) - {0}. 


(3) 


Equation (3) is recognlxad as a general eigenvalue problem of 
order M vitTi eigenvector (mode shape) and eigenvalue (cir- 
cular frequency squared). There are N values of u>^ and {^] so 
that the expansion of Equation (3) to include all solutions can 
be expressed as 


* r . * 

I I ‘**2 i 


|i>) 


2 

1 




“ [SI 





(4) 


u)2 

N 


or 


[A] [♦] - IS] l^]. 


(4a) 


Assume the structure Is grounded so that the stiffness matrix is 
nonsingular and its inverse exists. The inverse of the stiffness 
matrix is defined to be the flexibility matrix^ that is» 

[El = [SlT^ (5) 

Using this relationship in Equation (4a) and postmult ip lying by 
r -J gives 

[El [A1 [*1 - 1*1 (6) 

as Cha elganvalue problem to be solved in this paper. If the 
•tructure le not grounded, the stiffness matrix, [S], is singular, 
and Squacion (5) Is not applicable. The flexibility matrix, [E] , 
can still ha obtained by other techniques (either difeCtly or by 
operation on the stiffness matrix) so that Equation (6) is Still 
Valid. These techniques will not be discussed here. Jacobi's 
method cannot be used directly on Equation (6) because even 
though [£] and [A] are symmetric, their product In general is not. 
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tbs first operation in the solution of Equation (fi) is to 
dacoapose (raferenca Subroutine DCOHl) the mass i&atrlx into tri-^ 
angular factors* that is 

[A1 - [01^ tUl (7> 

vbmtm [U] is an upper triangular natrix. Next, we define 

[Tl - lUl [*J (8) 

fr<» which 

(♦1 - [or^ m. (9) 

Using these three -equations in Equation (6) and praoultiplylng by 
[Ul gives 

tDl I?1 - [T] r X J (10) 

where 

[D] - [U] [E] [U]^ (10a) 



[D] is referred to as the dynamic matrix in the subroutine* Be- 
cause the flexibility matrix [E] is symmetric, [D] is symmetric. 
Equation (10) is the eigenvalue problem from which the eigenvalue 
['X J and eigenvectors [T] are calculated by Jacobi's method 
(reference Subroutine ElGNl). The product of equation (9) to 
obtain the node shapes [^] is accomplished by using as the 

initial eigenvector in Subroutine EIGMl, thus eliminating a 
natrix multiplication. The i^^ circular frequency squared 
is obtained from Equation (10b) as 



( 11 ) 


Two orthogonality checks of the node shapes are calculated (on 
option) in the subroutine as follows. First, a summary of the 
T 

results of [4] [A] [<(] are printed. This calculation checks the 

T 

orthogonality of the mode shapes on the mass matrix (l.e. {(]>}^ [A] 
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and tha normalization ^i.e [A] » 1^. Next the results 

of [a] [K] [a] [0] are compared with (They should 

be equal*) This check results from premultiplying Equation (6) by 
T T 

[<&] [A] and assuming [^] [A] [^] - [I]* Only the upper half of 

[A] and [E] are used in the calculation of the mode shapes and 
frequencies but all of [A] and [£] are used in the orthogonality 
checks* Thus errors in the checks can result from nonsymmetric 
[Al or {El- 
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Subroutine MODEIX calculates the mode shapes and natural fre- 
quencies of a structure using its mass and stiffness matrices. 
This subroutine is a modification of Subroutine MODEl to allow a 
nonpositive definite mass matrix, ranove orthogonality checks, 
use 0 )^ convergence tolerance, o), f not calculated, and removal o'" 
making first element of each mode positive* The method of Jacobi 
is used in the solution of the general eigenvalue problem 




where the mass matrix is given by [A], the stiffness matrix by 
[S], and the size of the matrices by N* The mass and stiffness 
matrices must both be real and symmetric* The i,i element of the 

calculated diagonal matrix is the i^^ circular frequency 

squared* The corresponding mode shape is the i^^ column of [^1* 
Because Jacobi's method is used, all mode shapes and frequencies 
are calculated. The circular frequency squared, u)^, is output 
from the subroutine in increasing order of magnitude in {W2}. 

Th^ mode shapes are normalized (automatically due to Jacobi's 

method) so that [Aj = 1. 

DESCRIPTION OF TECHNIQUE 

Consider a discrete coordinate model of a structural system 
having N degrees of freedom. The undamped, free equations of 
motion for small vibrations about a position of stable equilibrium 
can be expressed in matrix notation as 


As before, [A] is the mass matrix and [S] is the stiffness matrix. 
(h(t)} is a column of discrete coordinate displacements. The 
time dependence in Eq [2] can be removed by the transformation 

{h(t)} ® e^^^ because the solutions of Equation (2) are 

assumed to be harmonic in time. Thus, for a nontrivial solution, 
Eq [2] becomes 

(-.u)^[A] + IS]) W = {0} . 

Equation [3] is recognized as a general eigenvalue problem of 
order N with eigenvector {4>} (mode shape) and eigenvalue oj'- 
(circular frequency squared). There are N values of and (4)} 
so that the expansion of Eq [3] to include all solutions can 
be expressed as 



lAl [41^ j {*}2 I 

I I 


or 


i 

I 
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I I 

IS] lu)^ I ; 

I I 


; {.y 


141 


[A] [♦] ^lu2g - is] [*]. 

Premultiplying by [A] ^ gives 


[4al 


lA]"^ [S] [$] = [®] [5] 

as the eigenvalue problem to be solved in this paper. Jacobi’s 
method cannot be used directly on Eq [5] because even though [A] 
and (SJ are symmetric, their product in general is not. 

The first operation in the solution of Eq [5] is to decompose 
the mass matrix into triangular factors, that is. 


[AJ = lU]’^ [U] [6] 

where [U] is an upper triangular matrix. Modified statements from 
Subroutines DCOMl are used, that is the square-root of absolute 
values are used. This allows for a nonpositive definite mass 

matrix. Using [A]~^ = ([U]^ [U])~^ = [U]"^ ((U]^)~^ and pre- 


mulciplying Eq 15] by [U} gives 

CtU]^)’^ [S] [$] = tu] [$] 17] 

Define 

[1] « lU] [®] (8] 

from which 

W ‘ [U]'^ KJ. [9] 

Using the fact that ([U]^) ^ = (tU]~^)^ we get 

[D] [*] = m [10] 

where 


ID] = au]"^)^ [S] [U]'^ 


tlOa] 
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and is referred to as the dyaamlc matrix in the subroutine. Be- 
cause the stiffness matrix [S] is symmetric, [D] is symmetric. 

Eq [10] is the eigenvalue problem from which the eigenvalues 

and eigenvectors [4>] are calculated by Jacobi's method (ref- 
erence Subroutine EIGNIA). The product of Eq [9] to obtain the 
mode shapes [4>] is accomplished by Using [U]~^ as the initial 
eigenvectors in Subroutine EIGNIA, thus eliminating a matrix 
multiplication. 



MULT 


Subroutine MULT calculates the product of two matrices. In 
matrix notation. 


where 


^^^NRAxNCB “ ^^^NRAxNRB ^^^NRBxNCB 


z 


ij 


NRB 

k^l 


( i = 1, NRa\ 
j = 1, ncb/ 


NRA is the number of rows of [A] and [Z). NRB is the number of 
rows of [B] and the number of columns of [A]. NCB is the number 
of columns of [B] and [ 3 ]. The number of columns of [A] must be 
equal to the number of rows of [Bl. 

Theorem : Multiplication of matrices is not commutati^^e in gen- 

eral. That is, 

[A][B] ^ [B][A] 


for any values of a.^ and b... 

ij ij 

Theorem : Multiplication of matrices is associative. That is, 

[A] ([B][C]) = ([A][B]) [C]. 

Theorem : Multiplication of matrices is distributive. That is, 

[A] ([B] + [C]) = [A][B] + [A][C], 


EXAMPLE 


Consider input of [A]j ^^2 ~ 


The reader can easily verify the output to be 


r 7. -8. 9.] 

[lO. 11. 12. J 


[Z]l,3 = [1. 2 


.] r 7. -8. 9.' 

[lO. 11. 12. 


= [27. 


lA. 33.]. 
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Subroutine MULTA calculates the product of two matrices. 
This subroutine is a modification of Subroutine MULT to allow 
larger matrix sizes by placing the answer [Z] in the same core 
locations as [A]. In matrix notation, 


where 


^^^NRAxNCB 


^^^NRAxNRB 


[B] 


NRBxNCB 


z 


ij 


NRB 

^ik ^kj 

k=l 


( i - 1, NRa\ 
j =« 1, NCB/ 


NRA is the number of rows of [A] and [Z]. NRB is the number of 
rows of [B] and the number of columns of [A]. NCB is the number 
of columns of [B] and [Z] * The number of columns of [A] must be 
equal to the number of rows of [B] • 

Theorem : Multiplication of matrices not commutative in gen- 

eral* That is, 

[A][BJ ^ [B][A] 


for any values of a, , and b, . . 

ij 

Theorem ; Multiplication of matrices is associative* That is, 
[A] (tBHCj) = ([A][B]) [C]. 

Theorem : Multiplication of matrices is distributive. That is, 

(A1 ([B] + [Cl) = [A][B] + [Al[C]. 


DESCRIPTION OF TECHNIQUE 


To accomplish matrix multiplication using only two matrix 
core spaces, an intermediate work space vector is used in the 
subroutine. The size of this work vector determines the limi- 
tation on the number of columns of [B] (i.e., NCB). The matrix 
multiplication is accomplished as follows. A single row of [A] 
is multiplied times the columns of [B] . These results are stored 
in the work vector until all the columns of [B] have been used. 
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The elements In the work vector then replace the row of [Aj used 
in the multiplication. This procedure is repeated for all the 
rows of [A]. Schematically, the procedure can je represented by 



where [Aq] is the original [A], [A^] » [Z] is the answer. 

EXAMPLE 

Consider input of [A], _ = [1. 2.] and [B] _ = j" 7. -8. 9.1 

|_10. 11. 12 .J 

The reader can easily verify the output to be 



[27. 14. 33.]. 
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Subroutine MULTB calculates the product of two matrices. 
This subroutine is a modification of Subroutine MULT to allow 
larger matrix sizes by placing the answer [Z] in the same core 
locations as [Bj. In matrix notation, 

^^^NRAxNCB " ^^^NRAxNRB ^®^NRBxNCB 


( i = 1, NRa\ 
j = 1. NCB^ 


where 


NRB 




ik kj 


k=l 


NRA is the number of rows of [A] and [Z]. NRB is the number of 
rows of [B] and the number of columns of [A] . NCB is the number 
of columns of [B] and [Z]. The number of columns of [A] must be 
equal to the number of rows of [B]. 

Theorem ; Multiplication of matrices is aot commutative in gen- 
eral. That is, 

[A][B] [B][A] 


for any values of a 


ij 


and b . . . 


Theorem ; Multiplication cf matrices is associative. That is, 
[A] ([B][C]) = ([A][B]) [C], 

Theorem : Multiplication of matrices is distributive. That is, 

[A] ([B] + [C]) = [A][B] + [A][C]. 


DESCRIPTION OF TECHNIQUE 


To accomplish matrix multiplication using only two matrix 
core spaces, an intermediate work space vector is used in the 
subroutine. The size of this work vector determines the limi- 
tation on the number of rows of [A] (i.e., NRA). The matrix 
multiplication is accomplished as follows. The rows of [A] are 
multiplied times a single column of [B]. These results are 
stored in the work vector until all the rows of [A] have been 
used. The elements in the work vector then replace the column 
of [B] used in the multiplication. This procedure is repeated 
for all the columns of [B]. Schematically, this procedure can 
be represented by 
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[A] 


1 

[Bq] 

1 1 

^work ) 
(vector/ 

[BlJ 

1 1 1 



1 

1 

1 1 
1 1 


I 1 1 

t 1 1 



1 

1 

1 1 
1 I 

■ ■ ^ {W} - 

1 {w}l I 

. 


1 

1 


J ' 1. 


where [Bq] is the original [B] , [BJ = [Z] is the answer . 

EXAMPLE 

Consider input of “ I-*-- “ T ^ -1 

[lO. 11. 12 .J 

The reader can easily verify the output to be 


IZ] 


1^3 


[ 1 . 



= [ 27 . 


14. 33.]. 
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Subroutine ONES generates a matrix with each element equal to 
one. That Is, 


In matrix notation, 

1 . 1 . ... 1 . 

1 . 1 . 

1 . ... 1 . 

where NR is the number of rows of [Z]; NC is the number of columns 
of (Zl. 
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Subroutine ONRBM calculates orthonormal ri 7 id body modes 
^'^^NxNRBM rigid body modes f f and mass matrix 

where N is the number of degrees of freedom and NRBM is 

.:he number of rigid body modes. The orthonormal rigid body modes 
orthogonal such that [A] {<{i = 0 (i?^j) where { 4 . ) is one 

rigid body mode (a column) of [4^] and i, j refer to che column 
number of [4>]. The orthonormal rigid body modes are normalized 
such that { 4 > 1 T[A){ 4 j}^ = 1 . 

DESCRIPTION OF TECHNIQUE 


Perform the triple matrix product 

= [B]. (1) 

Calculate the eigenvalues t > j and eigenvectors [E] of [B]. 

By definition. 


IE]^[B][E] = r X J 
or [B] = [E]-^f A HE]~1 

= [E]“^t Jt }[E]-} 

Substitution of Equation (3) into Equation (1), premultiplying y 
T 

[E] , and postmulLiplying by [E] yields 


( 2 ) 


(3) 


[E]^[^]^[A][^][E] = f /a Jt v'T j 


or 


J 


[E]^[^]^[Aj[.J.][E] 


w- 


[I]. 


The orthonormal rigid body modes are then defined as 


[4>] = [^][E] 


1 

/A 



PAGEHD 


Subroutine PAGEHD brings up a new page and prints a heading at 
the top of the page. This heading consists of: 

1) Run number; 

2} Date; 

3} Page number - Initialized as zero in subroutine START 
and incremented by one each time Subroutine PAGEHD is 
entered; 

4) User's name; 

5) Title Card 1; 

6) Title Card 2. 

Each of the above items was obtained in Subroutine START and trans- 
ferred by a COMIION block labeled LSTART. 
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Subroutine PLOTl plots points on a graph with linear x and y 
axes. The x, y coordinates of the points are supplied to the sub- 
routine in {XVEC} and [YMAT], respectively. Each column of [YMATJ 
is considered a separate curve with a maximum of three curves 
allowed. On option (IFCURV = 1), the points are connected by 
straight line segments. All curves will have the same y scale. 



PL0T2 


Subroutine PL0T2 plots points on a semilog or log-log graph* 
and connects the points with straight line segments « The input 
value of IPLOT determines the type of axes* IPLOT = 1 gives a 
log X axis and a linear y axis* IPLOT 2 gives a linear x axis 
and a log y axis* IPLOT = 3 gives both log x* y axes* The x* y 
coordinates of the points are supplied to the subroutine in (XVlXt 
and [YMAT], respectively. Each column of [YMAT] is considered a 
separate curve with a maximum of ten curves allowed* All curves 
will have the same y scale. 



PL0T3 


Subroucine PL0T3 plots points and connects the points with 

straight line segments to produce (on option) a perspective view 

or a stereo pair of a three^-dimensional object. The z co- 

t h t h 

ordinates of the i point of the object comprise the i row of 

the input matrix [CLOG). The points to be connected are specified 

in [MLOC] • Interpolation is done if a point falls outside of the 

cone of vision (usually about 60 deg) of the viewer* 



PLOTSS 


Subroutine PLOTSS selects the scale and calculates the top 
and bottom values of a 10-step linear scale from the imput values 
of the maximum and minimum values to be plotted. Each step (1/10 
of the total scale) will have a value of .5, 1, or 2. 
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Subroutine PSDl calculates power spectral density response and 
variance due to a power spectral density forcing function. The basic 
equation Is given in the frequency domain as 


( 


.2 



I q(") } * I F I u («a) 


where 


M is the mass matrix 
D is the damping matrix 
K is the stiffness matrix 
F is the force distribution 
u is the force variation 


( 1 ) 


The response | q(«) | could be obtained from 

jq(-)| » [m] + 1 - [d] + [k] )‘M F I «('-) (2) 

The disadvantage of this approach is that a complex matrix inversion 
must be performed for each « used. Because each matrix inversion re- 
quires much computer time, the approach of equation (2) is not used in 
this subroutine PSDl. Rather, the following approach used. 

The second order differential equation can be put into an equivalent 
first order differential equation given as 




- [mJ-‘Jd] ; - [m] [k] 



+ 

[m] |p| 

_{ q(t) 1 


H : i°] 


JqCOt 


]o} 


The response is expressed in terms of numerator and denominator containing 
products of first and second order polynomials permitting rapid evaluation 
at many • . Calculation of the complex roots uses considerable computer 
time but this technique is still faster than the technique of equation 
(2) when solutions at many • are sought. 


Define 



[- [m] ■ 

■‘W 

- [»] 


■ 




fj 


(3a) 
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lz(0| = 


! i(t) I 

I q(t) I 


, , W'MpI 

|b| 

|0| 

to give { Z(t) I = |^A*j { Z(t) I + I b I u(t) 


(3b) 


(3c) 

(4) 


The system characteristic matrix [,A*3 provides the basis for evaluating 
the resonant characteristics (natural frequencies) of the system des- 
cribed by matrices [m] , [dJ , and [k] . 

Taking the Laplace transform of equation (4) gives 

- (A*lj { Z(o)| = u I u(s) (5) 

Employ Cramer*s Rule to evaluate the pth element of |Z(s)| due to an 
input u, that is 


Zp(s) ^ Aug I Is - A* I 
u(s) |ls - A* I 


where Aug jls - A*|is accomplished by placing! b[ into column p of 

|ls > A*| . 

It is desired to evaluate both the numerator and denominator roots 
of equation ( 6 ). The Q-R algorithm (Reference J.G.F. Francis, ”The QR- 
Transformation - A Unitary Analogue to the LR-Transformation", The Com- 
puter Journal, Volume 4, October 1961, (Part 1) and Volume 5, January 
1962 (Part 2).) is a useful tool to extract the roots of the complex 
system and is used here. 

The denominator root extraction is straight forward in that we wish 
to find pj^, P 2 » P 3 , ... Pn from an expression of the form 


D(s) = det 




(7) 


such that 
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D(s) « (s - pj) (s - Pj) ... (s - pjj) 
N 

* TT ■ pp • 

i=l 


This evaluation is made by extracting the characteristic roots of the 
matrix M ■ In general) these roots will be complex because [A*] is 
not syometric. 

The numerator root extraction is more complicated in that the 
augmented matrix can very easily have a zero pivot element from | b | • 
ThuS) the null Pj^*s are eliminated from the expression giving the char- 
acceristic polynomial of order n less than N. It is a rather common 
occurence for the number of zeroes (order of the numerator N(s)) to be 
less than the number of poles (order of the denominator D(s)). The 
numerator expression N(s) can be written as 

n 

M, = kg TT 

i=l 


where the numerator root gain is 

n 

\ = (-!)“■" Tf Pi F* ■ ^""1 

i=l 


I is the identity matrix of size N with a null diagonal element. A^is 
the A* matrix with a negative column of jb[ in the pth column; and the 
Bode gain for the numerator is 

m 

s ■ <-*>” TT s^ where m^ n. (11) 

i=l 


Transfer function poles , zeros, and root gain can be converted to 
the standard Bode form for frequency resp nse by combining time con- 
stants (t) , damping (O and resonant frequencies («) as 
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vAiere the composite Bode gain is 


n 



Zi 


j 1, 


IT "j 

4=1 


( 12 ) 


(13) 


zi is a root of the numerator and, pj is a root of the denominator. 

The frequency response is then calculated by substituting jw for s and 
evaluating the transfer function expression (that is Zp(s)) at various 
M . In this subroutine the first and last « and the number of w wanted 
is specified by the user. 


The additional equations are obtained from 

^L(t)| = [aa] + [bb] <q(t)| + [cc] |q(t)| + | DD | 

or -- ^-<-2 [aa] + i<« [bb] + [cc] j + i {qi(“)| ) + {dd| 

which can be spile Into real and imaginary parts as 

jL^(“)| = [aa] + [cc] ^ |qj^(o.)} - [bb] |qj(«)| + { DD | 

[aa] + [cc] j |qj(“)| + [bb] {q^C")} 

Note that t q(*») } (real, imaginary) is the second half of | Z(<u) | , see 
equation (3b) . 
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Although the modal response Is calculated for each modal coordinate 
in turn for many frequencies, the need in the calculation of the addi- 
tional equations is for all of the modal coordinate responses at each 
frequency in turn. To avoid extensive tape reading, it was decided to 
place the modal responses into core for the additional equations cal- 
culation. Although this limits the number of frequencies that can be 
used, the loss in accuracy in calculating the variance should be minor. 

A constant step size between a specified minimum and maximum fre- 
quency is used. Because the area under the PSD-frequency curve (variance) 
is the final result, using a constant frequency step size gives a reasona- 
ble result even though some individual peaks may be missed. 

The PSD of the additional equations (at one station here) is given 
by 


L(w) = u(w) 


where H^(») = lJ(») + Lj(») 

and u(«) is the PSD-<»> variation of the forcing function | F ^ . The PSD-w 
variation is given in the input by TABF, TABW. 

The variance of the additional equation is 

"e 

f L(«) d» 

The square root of the variance (or mean square response) is the standard 
deviation. Printout (on option) in the subroutine includes: 

(1) - A* - 5 [a*] = 

(2) BGOL 5 I b I 

(3) R AR gives the real and Imaginary parts o£ the roots 

RARI gives the real and imaginary parts of the roots 
transpose. 

(4) RRED (printed in subroutine PSDl) 



-1 

1 

-1 1 

-M 

1' 

1 . 

1 ^ 

1 ^ 

I 

i 

1 

0 
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RRED(l) 

RRED(2) 

RRBD(3) 

RRED(4) 

RRED(5) 

RRED(6) 

RRED(7) 


is the nudber o£ numerator reals 

is the number of numerator complex pairs 

is the number of numerator zeroes 

is the number of denominator reals 

is the number of denominator complex pairs 

is the number of denominator zeroes 

is the Bode gain (Gg) > 


Based on the values of the first six quantities, come numerator 
time constants (r), damping (C) , and resonant frequencies (<•>), 
followed by denominator time constants, damping and resonant 
frequencies . 



PUNCHO 


Subroutine PUNCHO punches a matrix of octal numbers onto cards 
without round-off error. Octal representation of the matrix ele- 
ments is used because it gives an exact replica of the binary num- 
ber used by a digital computer. A decimal representation will 
not give an exact replica. The matrix on punched cards is to be 
used only as an emergency backup for the matrix written on a 
storage tape. 

This matrix on cards is compatible with the input form for 
Subroutine READO. A group of up to three consecutive elements 
from a row of "he matrix are punched on each card. If all of the 
elements of a group are zero, punching of this card is suppressed. 

The first card punched contains the matrix name (in card 
columns 1-6), the matrix row size (in card columns 7-10), and the 
matrix column size (card columns 11-15) . This is followed by the 
matrix data. On any card of the matrix data the first integer 
number (card columns 1-5) is the row number of the matrix elements 
on that card. The second integer number (card columns 6-10) is 
the column number of the matrix element in the first data field 
(card columns 14-25). The next group contains octal numbers (up 
to three numbers in card columns 14-25, 29-40, 44-55) that are 
the values of the matrix elements. This group of matrix elements 
is given in consecutive column order. The last card punched con- 
tains ten zeroes in card columns 1-10 to indicate the end of the 
matrix. 
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Subroutine RBTGi calculates a rigid body transformation 
matrix, C RBT J , such that 

HI - fsBiJ |q|REF 


where jq| is displacements at selected points on a structure, 
and IqjREF displacements at a reference point. The cartesian 
coorainates (x, y, z) of the selectco points are given by the 
input matrix UXYZ J . Each row of [ XYZ 3 has the x, y, z coor- 
dinates of one point. The x, y> z coordinates of the reference 
point are given byfXYZREPj. This reference point need not be 
one of the points given in fXYZj . A matrix of integers, £JDOFj 
is also input to this subroutine to define the degree of freedom 
number of the displacement of each point. This degree of free- 
dom number will define the row number of the degree of freedom 
intRBTj . Each of these degrees of freedom ar 2 assumed to be 
in the same direction as its corresponding reference degree of 
freedom. A negative value in t JDOF J indicates a relative (ver- 
sus absolute) degree of freedom (e.g., slosh, modal). The nega- 
tive value will cause that row of the resulting C RBTJ matrix to 
be zero. A row of integers, CJVECJ , is input to define the 
degrees of freedom associated with the reference point. Negative 
signs ontJVEC] enables change from assumed right hand system to 
one you wish to specify. Use of these matrices will be illus- 
trated in the example problem. 


DESCRIPTION OF TECHNIQUE 

th 

The translation (5^^, Sy, S^) at the i point is given in 
terms of the translation (5^^, Sy, S z) rotation 
of the reference point by the vector equation 




9 X V 
R i 


where Vi is the vector from the reference po?nt to the i*-“ 
point. All vectors are in terms of x, y, z components. 

The rotation at the i^^ point is equal to the rotation at 
the reference point, as shown oy the vector equation 
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A typical rigid body transformation for the point is 
size 6x6 and looks like 



where 

V = the X component of the vector distance from tk'^ 

reference point to point i, 

V = the y component of the vector distance from the 
^ reference point to point i, 

V = the z component of die vector distance from tha 

reference point to point i. 

The vector distance Is cal ulated by asing the i^^ row of the 
matrix fmj and the vector [XYZREFj , e,g., = XYZ (i, 1)- 

XYZREF(l). Using the Forma Subroutine REVADD, the 6x6 formed 
above is then **revadded** into the output matrix 
IVEC required for REVADD js formed from the i^^ row of the in- 
teger matrix DdofJ , and sign modified by the corresponding 
signs of the values in JVEC. The JVEC required for REVADD uses 
the integer vector CJVECJ . 

EXAMPLE 


A rigid body transformation will be calculated for the 
beam as shown in the sketch below that will express the rigid 
bedy motions of the 15 dof in terms of node paint 1 (the refer- 
ence point)* 
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Node 1 = 6 dof 

Node 2 = 3 dof 

Node 3 = 6 dof 

15 dof 


Input data: 

tX5fZ3 = A matrix of size NNODES (number c2 nodes) x 3 used 
to describe the coordinate locations of the node 
points. 


Node 1 
Node 2 
Node 3 


X 

0 . 

5. 

10 . 


y 2 . 

0 . 0 . 

7 . 6 . 

14. 12. 


tXYZREFj * A vector of size 3 used to describe the coordinate 
location of the reference point. The reference 
point is assumed to be Node 1. 


X y z 

Reference |^0. 0. oQ 

t JDOFj = An integer matrix of ' je NNODES x 6 used Lv des- 
cribe the degree-of-/^ dom table. In tl;is example 
there are no 0 deg; • :l of freedom at Node 2. 
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6y 

6z 

«x 

«y 


Node 1 

1 

2 

3 

4 

5 

6 

Node 2 

7 

8 

9 

0 

0 

0 

Node 3 

10 

11 

12 

13 

14 

15 


CJVEcJ An integer vector of size 6 used to define the 
degrees~of “freedom of the reference point. 

6x 5y 

Reference Q. 2 3 4 -5 £J 

Using the above input data, the rigid body transformation 
matrix, TRBXj , is calculated such that 

i^llSxl “^^’^^ISxb {‘^REfUxI* 

That is. 
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Subroutine RBTG2 calculates a rigid body transiormation 
matrix, CRBTJ , such that 

(,| - Cm3 1,1^ 


where 

|q| * displacements ( Sx» ®r> ®t* ^x» ^r» ^0 selected 
points on a structure 


and 


]ql_Ep * displacements (fix* *z» ®x» ^y» ^t a 
reference point. 

The following sketch is useful for identifying the variables 
involved. 
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of 

» y » 2 

coordina*:es of the reference point arc given by the input ma- 

trix [mREF] . 

The reference point need not be one of the points defined 
in Cxrt 3 . A matrix of integers, QdoiO , is input to the 
routine to define the degree^of- freedom number of the displace- 
ments of each point. This degree-of- freedom number will define 
the row number of the degree-of-freedom in . A negative 

value in indicates a relative (versus absolute) degree 

of freedom (e.g., slosh, modal). The negative value will cause 
that row of the resulting matrix to be zero. A row of 

integers, [JVEC^ , is input to define the degree-of-freedom 
associated with the reference point <i.e., the columns of 
Use of these matrices will be illustrated in the 
example problem. 

DESCRIPTION OF TECHNIQUE 

A cartesian rigid body transformation is formed first. 

The translations ( 8y, 8^) at the i^^ point are given in 

terms of the translations and rotations ( 8^, 8y, 82* ®x> ®y> 
of the reference point by the vector equation 

I = Sr + Sr ^ V. 


The cylindrical coordinates (x, 
points are given by the input matrix 

has the x, r, ^ coordinate of the 


CXIH^ ! 

1 th 


of the 
Tlic 
point , 


selcctc 
iil' row 
Tljc X 


where is the vector from the reference point to the i' 
point. The rotations at the i*^*^ point are equal to the 
tions at the reference point as shown by the vector equatic 1 



A typical cartesian rigid body transformation for the i^^^ 
point is size 6x6 and looks like 
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where 

= the X component of the vector distance from the 
reference point to point "i", 

» XRT(i,l) - XYZREF(l) 

V « the y component of the vector distance from the 
^ reference point to point "i", 

= XRI(i,2)*cos(XRl(i,3)) - ::xZREF(2) 

V ® the z component of the vector distance from the 

reference point to point "i", 

= XRT(i,2)*sin(XRT(i,3)) - XyZREF(3) 

Equation (1) may be written in matrix form as 


^%art^i " ^^^i {^artlllEF' 


( 2 ) 


It is now desired to transform jq^artli cartesian 

coordinates to cylindrical coordinates. This may be done by 
using a simple rotation transformation such that 
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r 




1 

«x 


1 1 
1 


Sx 

«r 


c s 1 

1 


«y 

*t 

i 

-s c * 
1 

1 

«2 

®x 


• 1 

\ 

»x 



c s 

1 





1 -S Cj 


e 

L tj 

1 





(3) 


vrfiere 


c = cos 4 ' = cos (XRT(i,3)) 
s = sin = sin (XRT(i,3)). 
Equation (3) may be expressed as 

Kylh = M hcartli 


Substituting equation (2) into equation (4) yields 

{'•cyltl ° I’catttREF 

■ HcartiREF <» 

Using the FORMA Subroutine REVADD, the resultant 6x6 
formed in equation (5) is “revadded** into the output matrix 
plBTj . The IVEC required for REVADD is formed from the i^^ 
row of the integer matrix £JD0F3 , and the JVEC required uses 
the integer vector C'^VEC] • 

EXAMPLE 


A rigid body transformation will be calculated fur the 
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ring, as shown in the sketch below, that will express the 
rigid body motion of the i8 DOF's in terms of the ircfercncc 
point. 

Node 1=6 DOF 
Node 2 = 3 DOF 

Node 3 = 6 DOF 

Node 4 = 3 DOF 

18 DOF 


Ring Geometry; 

X = 100. 

Diameter = 180. 

Node Points 90. Degrees 
Apart 

Input date: 

CxRTJ : A matrix of size NNODES (number of nodes)x3, used 

to describe coordinate locations of the node 
points. 

X r 

Node 1 FlOO. 90. 

Node 2 100. 90. 

Node 3 100. 90. 

Node 4 __100. 90. 

FxYZREfJ ; A scalar vector of size 3 used to describe the 
coordinate location of the reference point. 

X y z 

REFPT [30. 0. 30.3 

CjDOF 3 ; An integer matrix of size NNODESxb, used to 

describe the degrec-of-freedom table. In this 


90. 

180. 

270. 

360. 
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example, there ore no rotational clcgrees-of- 
freedom for node point 2 and no translational 
degrees-of* freedom for node point 4. 


Node 1 
Node 2 
Node 3 
Node 4 


1 

7 

10 

0 


8r 

2 

8 

11 

0 


DOF 


«t 

«x 

0 , 

<>t 

3 

4 

5 

6 

9 

0 

0 

0 

12 

13 

14 

15 

0 

16 

17 

18 


CJVECJ : An integer vector of size 6, used to define the 

degrees-of-freedora of the reference point. 

DOF 

^z ^2 

Reference Point £l 2 3 4 -5 oj 



the above input data, a rigid body 
, is calculated such that 


6x1 


transformation matrix 


1*^113x1 ~ t*^^18x6 HIrEF, 
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Subroutine READ reads a matrix of real numbers (a Fortran 
term for numbers with a decimal point) from either cards or tape 
into the computer. The matrix is then printed so that these in- 
put data are recorded with the answers of a run. A print sup- 
pression option is available for a matrix read from tape. On 
option, the matrix read from either cards or tape may be written 
on a tape (by Subroutine WTAPE) . 

The first data card read by Subroutine READ contains the in- 
formation to indicate whether cards or tape will be used. The 
Information entered on this card (and subsequent cards for card 
input) is given below. 


Card Data Input Form 


Required entries are denoted by an * symbol below. Any other 
entry is optional 


Card Format 

Columns Type (1) Entry 


First Card 1-6 


7-10 

11-15 

16-69 

72 


e 

o 


a 

o 


<u 

Cu 

<0 

0 ) 


u 


72 


73-78 


73-78 

73-76 


73-78 

79-80 


A ^Matrix Name. Will appear 

in printout. 

I ^Matrix Row Size. 

I *Matrix Column Size. 

A Any remarks to further iden- 

tify the input matrix. 

$. Only if the Write-Tape 
is to be initialized by Sub- 
routine INTAPE. The Write- 
Tape identification will be 
from card columns 73-78. 
or Anything other than $ is the 
Write-Tape is not to be ini- 
tialized. 

A The Write-Tape identification. 

(e.g., T1234). Use with $ in 
card column 72. 

or REWIND. The Write-Tape will 

be rewound before being used, 
or LIST. The Write-Tape will be 

listed by Subroutine LTAI^E 
after the matrix has been 
written on the Write-Tape, 
or Anything else will he Ignored. 

1 The Write-Tape Number, (e.g., 21). 

or Blank if the matrix Is not to 
be written on tape. 
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Card 

Format 


Columns 

Type (1) 

Entry 

Middle Cards 1-5 

I 

*Row Number of matrix elements 
on card. 

6-10 

I 

^Column Number of matrix ele- 
ment in first data field. 

11-27 

E 

*First data field with matrix 
elements. (2) 

28-44 

E 

^Second data field with matrix 
elements. (2) 

45-61 

E 

*Third data field with matrix 
elements. (2) 

62-78 

E 

^Fourth data field with matrix 
elements. (2) 

Last Card 1-10 

1 

*Ten zeroes. 


Nrte (1) Format Type A allows any keypunch symbol. 


Format Type I allows only integer numbers right justified 
In the field. Format Type E allows only real numbers 
(a Fortran term for numbers with a decimal point) any- 
where in the field. 

Note (2) Only nonzero elements need be entered. 

As an example of card input to Subroutine READ consider the 
following matrix: 



lAX'CJj,, - 

1. 

0. 

3. 

0. 

6. 

5. 




0. 

2. 

4. 

0. 

0. 

0. 




0. 

7. 

0. 

0. 

0. 

0. 

• 

This matrix Is 

also to be written 

on 

tape 

number 

21 that is to 

be initialized 

and identified 

as 

TA334. 

Figure 

1 demonstrates 


how this information could be written on a coding form to facili- 
tate keypunching to cards. 



FORMA 



Figure 1 Example of Card Input for Subroutine READ 




22-27 


or 


READ— 4/6 


th an ^ symbol below. Any 
card is used for each matrix 


Entry 

*Name of matrix to be read 
from the Re ad -Tape. 

Zero. The Read-Tape will 
move forward from its present 
position and search to the 
end of the tape. If the 
matrix is not found upon the 
first end-of-tape encounter, 
the tape will automatically 
rewind and make one more 
pass. If it Is not found on 
the second end-of-ta,e en- 
counter, an error message 
will be printed and the pro- 
gram will stop. 

Minus the location number of 
matrix on the Read-Tape. Tape 
will be positioned at the be- 
ginning of the location speci- 
fied and then continue as 
described above for a zero 
in column 10. 

*The Read-Tape Number, (e.g., 11). 
If positive, the matrix read 
will be printed in the output. 

If negative, the matrix read 
will not be printed in the 
output. 

*Run number of matrix to be 
read from the Read-Tape. 

REWIND. The Read-Tape will 
be rewound before being used. 
LIST. The Read-Tape will be 
listed by Subroutine LTAPE. 
Anything else will be concia- 
ered as part of the remarks 
described below. 
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Card Format 

Columns 'i'yP® (i) Entry 

A Any remarks to further iden- 

tify the input matrix. 

$. Only if the Write-Tape 
is to be initialized by Sub- 
routine INTAPE. The Write- 
Tape identification - .11 
from card columns 73-78. 
or Anything other than $ if the 
Write-Tape is not to be ini- 
tialized . 

A The Write-Tape identification. 

(e.g., T1234). Use with $ in 
card column 72. 

or REWIND. The Write-Tape will 

be rewound before being used, 
or LIST. The Write-Tape will 

be listed by Subroutine LTAPE 
after the matrix has been 
written on the Write-Tape, 
or Anything else will be ignored. 

I The Write-Tape Number, (e.g., 21). 

or Blank if the matrix is not 

to be written on tape. 

Note (1) Format Type A allows any keypunch symbol. 

Format Type I allows only integer numbers right justi- 
fied in the field. 

As examples of tape input to Subroutine Read consider: 

Example 1- A matrix naiaed AB2 with run number of RUN-46 is vO 
be read from tape number 11 info the computer and 
printed. This matrix is also to be written on tc.pe 
number 22 that is to be initialized and identified 
as T4321. 

Example 2. A matrix named XYZ4 with run number of TKD is on tape 
number 13 twice. The first rime is at location 29 
and the second time is at location 54. It is desired 
to read the second matrix. 

Figure 2 demonstrates how these two examples would be written 
on a coding form to facilitate keypunching to card' 


28-69 

r 72 



I 73-78 


V 1 

-76 

-78 
-80 





Figure 2 Examples o( Tape Input for Subroutine READ 
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Subroutine READIM reads a icatrix of integer numbers from 
either cards or tape into the computer* The matrix is then printed 
so that these input data are recorded with the answers of a run* 

A print suppression option is available for a matrix read from 
tape* On option^ the matrix read from either cards or tape may 
be written on a tape (by Subroutine Wi'APE) . 

The first data card read by Subroutine READIM contains the 
information to indicate whether cards or tape will be used. The 
information entered on this card (and subsequent cards for card 
input) is given below* 


Card Data Input Form 


Required entries are denoted by an * symbol below. Any other 
entry is optional. 


Card Format 

Columns Type (1) Entry 


First Card 1-6 


7-10 

11-15 

16-69 

72 


m 

c 

o 


72 


o. 

o 


73-78 


S’ 

V 


Vl 

:* 


73-78 

73-76 


73-78 

79-80 


A *Matrix Name. Will appear 

in printout . 

1 ^Matrix Row Size. 

I *Matrix Cciumn Size. 

A Any remarks to further iden- 

tify the input matrix. 

$. Only if the Write-Tape 
is to be initialized by Sub- 
routine I!n*APE. The Write- 
Tape identification will be 
from card columns 73-78. 
or Anything other than $ if the 
Write-Tape is not to be ini- 
tialized. 

A The Write-Tape identification, 

(e.g., TI234). Use with § 
card column 72. 

or REWIND. The Write-Tape will 

be rewound before being used, 
or LIST. The Write-Tape will 

be listed by Subroutine LTAPE 
after the matrix has been 
written on the Write-Tape, 
or Anything else will be ignored. 

I The Write-Tape Number, (e.g., 21). 

or Blank if the matrix is not 

to be written on tape. 
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Card 

Format 


Columns 

Type (1) 

Entry 

Middle Cards 1-5 

I 

*Row Number of matrix elements 
on card. 

6-10 

I 

*Column Number of matrix ele- 
ment in first data field. 

11-15 

I 

*First data field with matrix 
elements* (2) 

16-20 

etc 

I 

*Second data field with matrix 
elements. (2) 

76-80 

I 

^Fourteenth data field with 
matrix elements. (2) 

Last Card 1-10 

I 

*Ten zeroes. 


Note (1) Format Type A allows any keypunch symbol* 

Format Tyne I allows only integer numbers right iusti- 
fied in the field. 

Note (2) Only nonzero elements need be entered. 

As an example of card input to Subroutine READIM consider the 
following matrixi 


‘“*<=^3x6 


1 0 3 0 6 5 
0 2 4 0 0 0 
0 7 0 0 0 0 


This matrix is also to be written on tape number 21 that is to 
be initialized and identified as T4334. Figure 1 demonstrates 
how this information could be written on a coding form to facili- 
tate keypunching to cards* 



FORMA 
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Figure 1 Example of Card Input for Subroutine READ I M 



22-27 


or 
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th an * symbol below, ^ny 
card is used for each matrix 


Entry 

*Name of matrix to be read 
from the Read-Tape. 

Zero. The Read-Tape will 
move forward from its present 
position and search to the 
end of the tape. If the 
matrix is not found upon the 
first end-of-tape encounter, 
the cape will automatically 
rewind and make one more pass. 

If it is not found on the 
second end-of-tape encounter, 
an error message will be 
printed and the program will 
stop. 

Minus the location number of 
mat. ix on the Read-Tape. 

Tape will be positioned at 
the beginning of the location 
specified and then continue 
as described above for a zero 
in column 10. 

*The Reac-Tape Number, (e.g., 11). 
If positive, the matrix read 
will be printed in the output. 

If negative, the matrix read 
will not be printed in the 
output . 

*Run number of matrix to be 
read ^rcm the Read-Tape. 

REWIND. The Read-Tape will 
be rewound before being used. 
LIST* The Read-Tape will be 
listed by Subroutine LTAPE. 
Anything else will be con- 
sidered as part of the remarks 
described below. 
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Card Format 

Columns Type (1) Entry 

A Any remarks to further iden- 

tify the input matrix. 

$. Only if the Write-Tape 
is to be initialized by Sub- 
routine INTAPE. The Write- 
Tape identification will be 
from card columns 73-78. 
or Anything other than $ if the 
Write-Tape is not to be ini- 
tialized . 

A The Write-Tape identification. 

(e.g., T1234). Use with $ in 
card column 72. 

or REWIND. The Write-Tape will 

be rewound before being used, 
or LIST. The Write-Tape will 

be listed by Subroutine LTAPE 
after the matrix has been 
written on the Write-Tape, 
or Anything else will be ignored. 

I The Write-Tape Number, (e.g., 21). 

or Blank if the matrix is not 

to be written on tape. 

Note (1) Format Type A allows any keypunch symbol. 

Format Type I allows only integer numbers right justi- 
fied in the field. 

As examples of tape input to Subroutine READIM consider: 

Example 1. A matrix named AB2 with run number of RUN-46 is to 
be read from tape number 11 into the computer and 
printed. This matrix is also to be written on tape 
number 22 that is to be initialized and identified 
as T4321. 

Example 2. A matrix named XYZ4 with run number of TKD is on tape 
number 13 twice. The first time is at location 29 and 
the second time is at location 54. It is desired to 
read the second matrix. 

Figure 2 demonstrates how these two examples would be written 
on a coding form to facilitate keypunching to cards. 
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Figure 2 Examples of Tape Input for Subroutine READIM 
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Subroutine READO reads a matrix of octal numbers from cards 
into the computer. The matrix is then printed side by side in 
both octal and decimal so that these input data are recorded v/ith 
the answers of a run. 

The primary purpose of Subroutine READO is to read a matrix 
from punched cards without round off error. The cards are punched 
by Subroutine PONCHO. Octal representation of the matrix elements 
is used because it gives an exact replica of the binary number 
used by a digital computer. A decimal representation will not 
give an exact replica. The matrix on punched cards is to be used 
only as an emergency backup for the matrix written on a storage 
tape. 

Because of the emergency backup nature of input data to this 
Subroutine READO, only cards are read. No tape reading or writing 
options are available. 


The information entered on the data cards is given below. 
Required entries are denoted by an * symbol. Any other entry is 
optional. 



Card 

Columns 

Format 
Type (1) 

Entry 

First Card 

1-6 

A 

^Matrix Name. (Will appear 


7-10 

I 

in printout.) 
*Matrix Row Size. 


11-15 

I 

^Matrix Column Size. 


16-69 

A 

Any remarks to further iden- 

Middle Cards 1-5 

I 

tify the input matrix. 

*Row Number of matrix elements 


6-10 

I 

on card. 

^Column Number of matrix ele- 


14-25 

0 

ment in first data field. 
*First data field with matrix 


29-4l 

0 

elements. (2) 

^Second data field with matrix 


44-55 

0 

elements. (2) 

*Third data field with matrix 

Last Card 

1-10 

I 

elements. (2) 
*Ten zeroes. 


Note (1) Format Type A allows any keypunch symbol. 

Format Type I allows only integer numbers right justi- 
fied in the field. Format Type 0 allows only octal 
numbers. 

Note (2) Only nonzero elements need be entered. 
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No examples of input are given because data would not be key- 
punched for input to Subroutine READO but rather obtained from 
Subroutine FUNCHO. 
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Subroutine REVADD rearranges (re 'ises; the rows and columns of 
a matrix [A]* multiplies [A] by a scalar alpha ^ and adds the re- 
sult to a previously defined matrix [Z]. The revision of [A] is 
specified by two vectors. The first vector |lVEc[ gives the new 
row location of each row of [A] in [Z). The second vector |jVEC| 
gives the new column location of each column of [A] in [Z]. The 
[Z] matrix must be defined before the use of this subroutine. For 
instance, if [Z] is to be originally defined as all zeros. Subrou- 
tine ZERO could be used. The REVADD operation can be thought of 
in subscript notation as 




(la) + a 


/I = 1, NRAv 
U = 1, KCA) 


where 


k = IVEC(l), 
i = JVEC(j). 

NRA is the number of rows of [Aj , and NCA is the number of columns 
of [A] . 

Values in |lVEC| and |jVEc| may be positive, negative or zero. 
A negative value changes the sign of the corresponding row or 
column of [A] in [ZJ. A zero value omits the corresponding row or 
column of [A] from [Z]^ The values are integer numbers. 

This subroutine may be called repeatedly to form [Z] from the 
revision/addition of several [A] matrices. 

An important use of Subroutine REVADD is to revise and add the 
stiffness matrix of a structural component to the stiffness matrix 
of the complete structure to account for the difference in coor- 
dinate systems. 

EXAMPLES 


The first example to illustrate the REVADD operation is as 
follows: 

Matrix [Z] has been previously defined as 


U) 


4x5 


1. 

0. 

0. 

0. 

0. 

0. 

2. 

0. 

0. 

0. 

0. 

0. 

3. 

4. 

5. 

0. 

0. 

0. 

0. 

6 . 
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Let [A] be defined as 

2 . 

4 , 

6 . 

The first row of [A] is to be added to the third row of [Z], the 
second row of [A] is to be omitted from [Z], the third row of [A] 
is to be added to the first row of [Z], with the sign of each 
element reversed « 



Thus {iVECf 


3x1 


3 

0 

-1 


The first column of [A] is to be added to the second column of [Z] 
with the sign of each element reversed, the second column of [A] 
is to be added to the fifth column of [Z]« 


Thus |JVEC| 2^1 


-2 

5 


Then, assuming 

u =» 

1.0 

and placing jiVEC 

[ and 

jjVECf adj 

acent 

to 


[A] to aid in visualizing the 

revision of [A] 

, we 

have 





1. 

0. 

0. 

0. 

0. 



1-^ 

i 






0. 

2. 

0. 

0. 

0. 

+ 1.0 

3 

1. 

2. 






0. 

0. 

3. 

4. 

5. 


0 

3. 

4. 






_0. 

0. 

0. 

0. 

6._ 


-1 

5. 

6. 






'l. 

0. 

0. 

0. 

0 " 









a 

0. 

2. 

0. 

0. 

0, 

+ 1.0 

3 

0. 

-1 


0, 

0. 

2. 


0. 

0. 

3. 

4. 

5. 


0 

0. 

-3 


0. 

0. 

4. 


0. 

0. 

0. 

0. 

6. 


-1 

0. 

-5 


0. 

0. 

6. 
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A second example of the use of this subroutine demonstrates 
the coordinate transformation concept which is a very important 
application of REVADD. It is desired to transform a stiffness 
matrix from an original coordinate system (subscript 1) to a 
final coordinate system (subscript 2) as shown in the sketch below. 





For this example, the stiffness matrix in the original coordinate 
system is assumed to be 
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From inspection of the coordinate 8ys*-em axes in the above, sketch, 



To obtain the vector on the right hand side of the equation: 

xi of the original coordinate system is to be the second 
row of the final coordinate system with the sign reversed. 

Yl of the original coordinate system is to be the third 
row of t.ie fin«l coordinate system. 


zj of the original coordinate system is to be the first 
row of the final coordinate system with the sign reversed. 

The |lVEC| to accomplish this change is 


3x1 


Apply this |lVEc| as the |lVEc| and |jVEc[ to the original stiff- 
ni-s:; matrix to obtain the final stiffness matrix. This application 
i*- analogous to the change in coordinates r.^Je in the triple matrix 
product procedure. 


t- ■ -] 





The above REV ADD procedure ir.ay be compared with the more conven- 
tional triple matrix produce procedure below. 
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The equation £or the coordinate transforniation v-ould be 



The strain energy expression with the stiffness matrix is 
T 

fwiveu by U * ^2 I'K] |qj. or for coordinate system 1, 



Substituting Equation (1) into (2) gives 


( 1 ) 



thus the inner triple matrix product gives the stiffness matrix 
in coordinate system 2. i,e.. 



which is the same result as that obtained from the REVADD procedure. 

The advantages of using the REVADD procedur-*: over the triple 
matrix product procedure are: 

1) Less computer time; 

2) Less computet core ?s used; 

3) Usually easier to code the jlVEcJ (thus |jVEC|) 
than to code the transformation matrix. 



ROUMLT 


Subroutine ROUMLT evaluates a special matrix operation by 
multiplying each row of a matrix [B] by a scalar. That is. 


U1 


MRxNC 


w’ “ire 


aj {bj} 


IxNC 

T 


NR 


I**nrI 


IxNC 


a. . 
ij 


i. b.. 
1 ij 


/i » 1. nr\ 

\i * 1. Ncy 


Each scalar a^ is an element of the 


|b^| denotes row i of [B]. 

input vector {AVEC}. NR is the number of tows in [B] and [Z] and 
the size of {AVEC}. NC is the number of columns in [B] and [Z] . 
The nuober of elements of {AVEC} must be equal to the number of 
rows of [B]. 


EXAMPLE 


Consider the input of 


(AVEC)j^j 

The output will be 
IZJ 




and (B] 


2x3 


[ 7. 2. 

4. 5. 


2x3 



1 . 17. 

2. 

-3 

L- 

J. [4. 

5. 

6 

r 14. 

4. 

-6.} 

[-12. 

-15. 

-18.J • 




Subroutine RTAPE reads a sel».:ted matrix from tape (disk) into 
the computer core. The matrix tp be selected is identified by the 
desired run number and matrix i\<^e. This procedure is accomplished 
by searching the matrix headings (see Subroutine VITAPE vnriteup) 
until a match with the desired run number and matrix name is ob- 
tained and then reading the matrix elements from tape (disk) into 
Che cixnputer core. The search starts from the current position 
(does not rewind) of the tape (disk) and proceeds to the EOT (end 
of Cape defined in Subroutine WTAPE writeup). If the desired 
matrix was not found upon reaching the EOT, a rewind is performed 
and one more search "o the EOT is made. If the desired matrix is 
again not found, (1) an error message is printed, (2) a listing 
of Che matrix headings is printed (see Subroutine LTAPE writeup) , 
and (3) transfer is made to Subroutine ZZBOMB where the program 
is terminated. 



SIGMA 


Subroutine SIGMA generates a square matrix with elements on 
and below the diagonal equal to one and above the diagonal equal 


Co zero. That is. 



'ij ■ 

a 

iJ) 

hi ■ 

(i 

< j) 


In matrix notation, 


1 . 1 . 


• ♦ 

1 . . . . 1 

where N is the size of [ZJ « 

This subroutine is useful in calculating the sum of all quan- 
tities ox a given collection* The name of the subroutine comes 
from the algebraic notation for this process, that is, Z* 

As an example of the use of this subroutine, consider the fol- 
lowing beam with leads p^ at discrete points on the beam: 



The shear to the right of each point (i) on the beam is 

or in matrix notation, 

Vj 1. 0. 0* 0. ”Pi" 

V2 1* 1* 0. 0* P2 

V3 “ 1. 1. 1. 0. P3 

V14 1* 1. 1. 1. P4 
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Subroutine SMEQl gives the solution of linear simultaneous 
algebraic equations. These equations are expressed in matrix 
form as 


(Al 


NxM 


{Z}„ , = - 

Nxl Nxl 


where [Aj and {B} are known and {Z} is to be calculated. N is 
the size of the system. The method, with which Gauss's name is 
associated, consists or a chain of successive eliminations by 
which the original system is transformed into a system with an 
upper triangular matrix whose solution is easily calculated. A 
modification of the standard technique is used in the subroutine. 
That is, the largest pivotal divisor is searched for and used. 

The rows are interchanged when necessary to accomplish this. This 
technique will give the most accurate results since division by 
small numbers will be avoided. 


DESCRIPTION OF TECHNIQUE 

Given a system of "N" linear algebraic equations, 

ail ^1 *12 ^2 *1K *N “ *^1 

*21 *1 *22 *2 *2N *N * ^2 

♦ 

* 

*N1 *1 *N2 ^2 ‘ *NN ^N * *’n’ 

The procedure used to obtain {Z} is as follows: 

1) The largest coefficient of Zi is found. The equation in 
which this occurs is interchanged with the first equation; 

2) Divide the new first equation by the new 

3) Alter Equations 2 through N by subtracting from the i^^ 
equation the value Sn times the first equation from step 

2). This is done for i = 2, 3, . . . N. 
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After these three steps* a new but equivalent system of equations 
resiits: 


*1 *12.1 *2 *13.1 *3 ®iN.l ” **1.1 

*22.1 *2 *23.1 *3 *2N,1 *N “ ^2.1 

*N2.l *2 *N3.1 *3 + • • • ■*■ *nn.1 *N * \.l 

where the subscript (»1) designates the first elimination. The 
variable has thus been eliminated from Equations 2 through N. 

Steps 1) , 2) , and 3) above are repeated on Equations 2 through 
N such that z^ is also eliminated. This process is continued un- 
til an equivalent system of equations is of triangular form as 
shown below. This process is referred to as the forward solution. 


— 



“ — 


• *■ 

^ *12.1 *13.1 

*1N.1 


"l 


^.1 

1 

*23.2 

• 

A 

*2N.2 

1 


*2 

• 

• 

* 

2.. 

= 

*’2.2 

• 

e 

• 



1 

N 

. j 


N.N 


The back solution for {Z) from the above system of equations 
gives no problem since z^ is immediately given as b^^ Next 

z^ ^ is easily found from the (N-1) Equation and so on until all 

z^ are known. 
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EXAh?LE 


The following data are given as input. 


^^^3x3 

2. 

5. 

3. 

6. 

4. 

7. 

and = 

2. 

-3. 


8. 

9. 

9. 


4. 


Thus the equations to be solved are 


2. 3. 4. 




2. 

5. 6. 7. 


^2 

II 

-3. 

8. 9. 9. 


/3. 


I 

• 

L. 


Using the technique described previously, the reader can verify 
the output to be 




■-19.00' 

^2 

= 

29.33 

/3. 


-12.00 
. - 


REFERENCE 

Faddeeva, V. N.: Computational Methods of Linear Algebra. 

Dover Publications Inc*, New York, 1959. 
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Subroutine SRED2 operates on the stiffness matrix [A] to form a 
reduced stiffness matrix [R] and/or the reducing transformation 
[T]. The relation between the stiffness matrix [A], displacements 
(X), and applied forces {B} may be expressed in matrix form as 

[A] {X} - {B} tH 

The reduction method assumes Eq [1] to be partitioned as 

[All] 

1^21 ^^ 22 ^ 

where {X^} are the displacements to be reduced out and {X2} 

are the displacements to be retained* The applied forces acting 
on the coordinates to be reduced are assumed to be aero» such that 

=» {0} [3] 

Substituting Eq [3] into Eq [2J and expanding the upper partition, 
will yield the reduced displacements in terms of the retained 
displacements as 

{x^} - - ^*2^ 

Expanding the lower partitions of Eq [2] and substituting Eq [A] 
will yield the reduced stiffness matrix as 



[R] = {B^} [5] 

where [R] is the reduced stiffness matrix and is expressed as 

IR] - IA22I - lA 2 jl lA^^r^ [A^ 2 ^ 

The reducing transformation [T] may be expressed using Eq [4] as 


V2U 


[T] {X2> 


[7] 



where 
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m • 

Also 




m 


[R] = [T]^ [A] [T] 


181 


19 ] 


DESCRIPTION OF TECHNIQUE 

This subroutine uses Gauss reduction partially completed to form 
matrix [R] and [T] from stiffness matrix [A]. As an example of 
the method, consider three simultaneous equations of the following 
form 


*11 ^1 *12 *2 *13 *3 “ ^1 
*21 ’‘l *22 ^^2 *23 *3 “ ^2 
*31 ^1 *32 ^^2 *33 ^3 “ ^3 


These equations may be written in matrix form as 


*11 *12 *13 


*21 *22 *23 


*31 *32 *33 



Solve the first Equation for as 


*12 

r ~ ^^2 


“13 ^ h 

*11 ^ *11 


un 


1121 


Substituting Eq [12] into ' e second and third equations in Eq [11] 
and divide the first by a,, results in 


1 

A 

®12 

* ~ 
*13 


''Xi”> 


ro 

0 

A 

*22 

* 

*23 

< 

^2 

r< 

^2 

0 

it 

*32 

* 

*33 






[ 13 ] 
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where 



* 

h2 


h2 “ 

hi 


* 

h3 


h3 “ 

hi 


* 

^22 “ 

*22 ~ 

*12 *21 
*11 

* 

A S 

23 

*23 

*13 *21 

hi 

* 

^32 “ 

*32 ■ 

*12 *31 
*11 

* 

h3 “ 

*33 

*13 *31 

hi 


H * '’1 
b = 


bz = b2 - 


b3 = 63 - 


^21 

hi 

h hi 

hi 


Solve the second equation for x- which will yield 

* ^ 


V a: — 
*2 


®23 


*22 


^3-^ IT 

h2 


Substitute Eq (23] into the third equation in Eq [13] and divide 
the second equation by a^z* This results in 


* 

*12 


a 


1 

0 


13 

** 

*23 

*33 





r* 


h 


h 

< 

h 

Jl^ 

** 

**2 


X 

V. -V 


b** 

J 


114] 

115] 
[16] 

[17] 

[18] 

[19] 

[ 20 ] 

[ 21 ] 

[ 22 ] 

[23] 


[ 24 ] 



where 
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* 


itik 

*23 


23 * 

it 

*22 

it it 

** 

it 

*23 *32 

CJ 

It 

*33 " 

it 

A 


*22 


** ^2 


n 

CM 

* 

*22 


itit 

it 

k — 

it it 

^2 *32 

3 

^3 

it 


3 a* 
*22 


[25] 

[26] 

[27] 

[28] 


The reduced stiffness matrix has been formed and is contained as the 
** 

a_- element. This can be shown if in Eq [2] we let 


{X2} = {x^} 



[A21] » [a^j^ 


^^22^ “ *33 




[31] 

[32] 

[33] 

[34] 

[35] 


[ 36 ] 
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Substituting Eq [31] through [34] into Eq [6] results in a reduced 
stiffness matrix of the form 


[R] = 


^33^ 


^21 ^23 ^31 ^12 ^13 ^32 ~ ^13 ^ 2 .! ^31 ~ ^11 ^23 ^32 

®11 ®22 ■ ®12 ®21 


[37] 


Equation [37] is identical to the result obtained by expanding 
Eq [26]. Thus, Gauss reduction partially completed yields the 
reduced stiffness matrix. 

The reducing transformation may also be obtained using Gauss re- 
duction if additional operations are performed. From Eq [24], 
solve the second equation for X 2 


icic 

*2 “ “ *23 *3. ^2 

Substitute Eq [38] into the first equation in Eq [24], which will 
yield 


kkit 


r 



1 0 aj3 





kit 




kk 

0 1 a23 

< 

^2 

1 '^ 

A* 




** 1 

0 0 a,. 





33 

— 


3 

J 


L 3 j 


where 

*** it it irk 

*13 “ *13 ■ *12 *23 


*** * * ** 
'’l “ ^ ' *12 ^^2 


[38] 


[39] 


[40] 

[41] 


Inspection of Eq [39] shows that we have formed a unity matrix 
partition where the row-columns were reduced. The a** and the 

kit 

022 elements contain the necessary information to form the reduc- 
ing transformation. To show this, substitute Eq [31] and [32] into 
Eq [8] to yield 
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*21 ^23 ' ^13 ^22 
*11 * 22 ." *21 *12 


[T] 


*12 *13 " *11 *23 
*11 *22 " *21 *12 


1 


[42] 


The second and third rows of Eq [42] are equal to the negative of 

icitit irk 

elements and a^^ In Eq [39]. Thus, Gauss reduction also 

yields the reducing transformation. 



START 


Subroutine START performs the following operations: 

1) Reads Input Card 1 for the run numbei (any keypunch 
symbol in card columns 1 thru 6) and 'he user's name 
(any keypunch symbol in card columns 11 thru 28). 

If the run number is equal to STOP (Card columns 1 
thru 4), the run is terminated. 

If the run number is not equal to STOP, the run con- 
tinues in Subroutine START as follows. 

2) Reads Input Card 2 for Title Card 1. Any keypunch 
symbols may he used in Card columns 1 thru 72. 

3) Reads Input Card 3 for Title Card 2. Any keypunch 
symbols may be used in Card columns 1 thru 72. 

4) Initializes page number as zero for use in Subroutine 
PAGEHD. 

5) Interrogates computer for the date. 

Run number, date, page number, user’s name, Title Card 1, and 
Title Card 2 are transferred by a COMMON block labeled LSTART for 
use in other subroutines PAGEHD, PLOTx, PL0T2, PL0T3, and WTAPE, 

Subroutine START is used to start each computer run in the 
FORMA system and will normally be the first subroutine called in 
a computer program. As an example, pertinent statements from a 
program using S - RT could be: 

1 CALL START 


GO TO 1 
END 
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Subroutine STIFl takes distributed longitudinal stiffness of 
a rod and replaces it with a *'free-free'* stiffness matrix. Lon- 
gitudinal stiffness will be discussed here but the results are 
also applicable to the torsional stiffness case. Tlie elements of 
the stiffner.s matrix are representative stiffness values of the 
rod at selected points ^n the rod. These elements are calculated 
by assuming constant axial force between pairs of the selected 
points • 

The x-stations of the selected points (panel points) are given 
in {PP}. These x-stations must be in increasing order. 

The distributed stiffness » A£(x) » is assumed to be piecewise 
linear and is represented by straight line segments as shown in 
Figure 1. 

A£(x) 



The x-stations of the first and last points for distributed 
stiffness must coincide with the first and last panel point x- 
stations, respectively. All other x-stations of the end points 
for the line segments giving the distributed stiffness are inde- 
pendent of the panel point x-statlons. The line segments may or 
may not be joined; however, there must not be any x voids or over- 
laps. The distributed stiffness is defined in [DAEJ. Each row 
of [DAE] represents one nonvertical line segment. The form of 
each row of [DAE] is [x^ X 2 AE^ AE 2 I where , A£| give the first 
end point and X 2 » AE 2 give the second end point of a 3ine segment* 
The X 2 of row 1 of [DAE] must be equal Co of row 2, etc. 

The calculated representative stiff ▼'ess values at the selected 
panel points are placed in a stiffness matrix given by [Z] which 
is symmetrical and Cri-dlagonal. That is. 
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*1.2 

*2,2 *2,3 
*3,2 *3,3 *3,4 

n-1 , n-2 n-1 , n-1 n-1 , n 

z 1 z 
n»n-l n,n 

vhere z • z . and n is the number of panel points. The gen* 

eralized coordinates associated with [Z] are axial deflections at 
each of the panel points • 

As mentioned before» the results of the longitudinal stiffness 
case considered in this paper are also applicable to the torsional 
stiffness case. The following table gives the relationship of 
variables 


Longitudinal 

Torsional 

X 

X 

6 


AE 

GJ 

P (axial force) 

T (torque) 


[Z1 


Qxn 


‘ 1.1 

‘ 2,1 


DESCRIPTION OF TECHNIQUE 

The replacement of distributed longitudinal stiffness of a 
rod by a stiffness matrix is obtained using a strain energy ap- 
proach as follows. Consider a portion of a rod that is loaded 
with an axial force at panel point k and restrained at panel 
point k + 1 as shown in Figure 2. The region between panel points 
k and k + 1 is referred to as bay k. 



Figure 2 Loading on Bay k 
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Tha strain energy In bay k is defined by 


U. = 


1 


^+1 


P(x)- 


A(x) E(x) 


dx« 


( 1 ) 


where 


P(x) is the axial force* 

A(x) is the cross-sectional area of the rod, 

E(x) is the modulus of elasticity of the rod material, and 
X is the longitudinal axis of the beam. 

To integrate Equation (1) , the axial force in bay k is assumed 
constant and equal to the axial force at panel point k. That is* 

P(x) = P^. (2) 

Also, the product A(x) £(x) is assumed to vary linearly as shown 
in Figure 3. 



Figure 3 Stiffness Distribution 


The equation for a straight line segment as shown in Figure 

3 is 


AE(x) = AEi + (x - xi) (AE 2 " AEi)/(x 2 “ xi). (3) 

Substituting Equations (2) and (3) into Equation (1) gives the 
strain energy of the axial stiffness represented by one line seg- 
ment i In bay k as 
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( 4 ) 


The subscripts p and q have been introduced to handle the possi- 
bility of a line segment extending past the bay limits* Thus» x 

P 

is the greater of x. or x. and x is the lesser of x^ or x, . ^ 

® 1 K q 2 k+1 . 

Similarly, AE^ is either A£^ or AEj^, and A£^ is either A £2 or 
AEj^^l* The Integration is continued for the line segment in ad- 
jacent bays, if necessary, until the entire line segment has been 
used* Performing the integration of Equation (4) yields 

U. . = f 
i.k k k 


where 


1 ^ 

p 


“pI/I*, - ‘p) • 


For constant stiffness, i.e., AE^ = Equation (5a) is of in- 

definite form* For this case* integration of Equation (4) yields 


^k = 


(x - X ) 

^ q - pJ 

AE 


(5) 

(5a) 

(5b) 


(5c) 


Equations (5) are evaluated for every line segment of distrib- 
uted stiffness. Then, application of Castigliano's Theorem gives 
the axial deflection of panel point k relative to panel point 
k + 1 as 


from which 


A6. 


^"k 

— - = P f 
3P, k k 
k 



( 6 ) 
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The restr aint at panel point k + 1 is removed by application of 
the era Information 



11 -11 



(7) 


Sul stitut Lon of Equations (6) and (7) into Equation O) (evalu- 
ated foz.' all line segments 1 in bay k) gives the strain energy in 
bay k as 


^k+l] 

where 

\.k “ \+l,k+l “ V^k’ 
*k,k+l “ "V^k* 



k.k 

(sym) 


k,k+l 


k+l,k+l 



' 6 . 


k 


"^k+l 


( 8 ) 


(8a) 

(8b) 


The ke:^neI matrix in the triple matrix product of Equation (8) is 
the stiffness matri> »:hat represents the longitudinal stiffness 
in bay k. 

The stiffness matrix for the entire rod is obtained by evalu- 
ating Equations (8), (5a), and (5b) for every bay. Each resulting 
2x2 bay stiffness matrix is added to previous 2x2 bay stiffness 
matrices at like panel points to form the stiffness matrix (free- 
free) for tha entire rod. 


SPECIAL CASE 


For constai, stiffness (A£) extending from x 
the stiffness matrix is 


P 


to x 

q 


in bay k. 
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Subroutine STIF 2 takes distributed bending stiffness and (on 
option) dlstrlbut*.d shear stiffness of a beam and replaces them 
with a "free-free" stiffness matrix. The elements of the stiff- 
ness matrix are representative stiffness values of the beam at 
selected points on the beam. These elements are calculated by 
assuming constant shear and linearly varying bending moment be- 
tween pairs of the selected points. 

The x-stations of the selected points (panel points) are given 
in {PP}. These x-statlons must be In Increasing order. 

The distributed bending stiffness, EI(x), is assumed to be 
piecewise linear and Is represented by straight line segments as 
shown in Figure 1 . The x-statlons of the first and last points 
for distributed stiffness must coincide with the first and last 
panel point x-statlons, respectively. All other x-statlons of 
the end points for the line segments giving the distributed stiff- 
ness are Independent of the panel point x-statlons. The line 
segments may or may not be joined; however, there must not be any 
X voids or overlaps. The distributed bending stiffness is defined 
in [DEI]. Each row of [DEI] represents one nonvectlcal line seg- 
ment. The form of each row of [DEI] is [xj X2 EIj EI2] where xj , 
EI^ give the first end point and X2, £12 give tiie second end point 
of a line segment. The X2 of row 1 of [DEI] must be equal to X]^ 
of row 2 , etc. 


EI(x) 



Figure 1 Distributed Bending Stiffness 

The distributed shear stiffness, KAG(x), is also assumed to 
be represented by straight line segments. All statements used 
above for distributed bending stiffness are applicable to distrib- 
uted shear stiffness. The distributed shear stiffness is de- 
fined in [DKAG]. The form of each row of [DKAG] is [xj X2 KAGi 
KAG2]. The end point stations x^ and X2 for the shear stiffness 
line segments are independent from the end point stations xi and 
X2 for the bending stiffness line segments. 
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The calculated representative stiffness values at the selected 
panel points are placed in a stiffness matrix given by [Z]. The 
form of [Z] is 


[Z1 


2nx2n 




1 

1 

• 

_ j 

^.1 


iv.] 

i M 


where n is the number of panel points* Matrix [Z] is symmetric, 
i*e., *» ^ji’ partition fonu of [Z] results from the two 

generalized coordinates (lateral translation 6 and rotation 0) at 
each panel point arranged with all translation coordinates first 
followed by all rotation coordinates* Each partition of [Z] is 
square and tri-diagonal, for example. 


"" 1,1 ^ 1,2 

^ 2,1 ^ 2,2 ^ 2,3 



nxn 


^ 3,2 ^ 3,3 ^ 3,4 


^n-l,n-2 ^n-l,n-l ^n-l,n 
^n,n-l ^n,n 

The sign convention used in this paper to obtain the stiff- 
ness matrix [Z] is shown in Figure 2. 


4 


6 


e 




Figure 2 Sign Convention 
for This Paper 
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Th« beam stiffness matrix obtained In this paper is applicable 
for either the pitch or yaw plane with a change of variables and 
possible ?lgn changes. For the axis system shown In Figure 3, the 
followirg variables would be used. 


This 
Paper 
(6, X, 0) 

Axis System of Figure 3 

Pitch 

("■ *■ M 

Yaw 

(y. X. 0,) 

X 

X 

X 

6 

6 

6 


z 

y 

e 

0 

-0 


y 

z 



Figure 3 Beam Axis System (Right-Hand) 

Note that in order to have a right-hand system such as that shown 
in Figure 3, the sign of the and partitions of [Z] 

.rom this subroutine would have to be changed for the yaw plane. 
DESCRIPTION OF TECHNIQUE 

The replacement of distributed bending and shear stiffness 
of a beam by a stiffness matrix Is obtained using a strain energy 
approach as follows. Consider a portion of a beam that Is loaded 
with a shear and moment at panel point k and restrained at panel 
point k + 1 as shown In Figure 4. The region between panel points 
k and k + 1 Is referred to as bay k. 
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Figure 4 Loading on Bay k 


The strain energy in bay k is defined by 



where 

M(x) is Che bending moment, 

V(x) is the shear, 

E(x) is the bending modulus of elasticity of Che material, 

l(x) is the cross-sectional moment of inertia about Che 
beam's neutral axis, 

K is the shape factor (e.g., K = 1 for a solid circular 
cylinder, K ■ 0.5 for a thin walled circular cylin- 
der) , 

A(x) is Che cross-sectional area, 

G(x) is Che shear modulus of elasticity of the material, 
and 

X is the undeformed longitudinal axis of the beam. 

To integrate Equation (1), the following assumptions are made. 
First, the shear in bay k is assumed aonetant and equal to the 
shear force at panel point k, that is, 


V(x) - V^. 


( 2 ) 
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Second t the bending moment in bay k is assumed to vary lineax'ly ^ 
that ISt 

M(x) = - Xj^j. (3) 

Third, the products E(x) l(x) and KA(x) G(x) are assumed to vary 
linearly as shown in Figure 5. The equation for a straight line 
segment as shown in Figure S is 

EI(x) = Ell + (x - xi) (EI 2 - EIi)/(x 2 - xi). (4) 


EI(x) 



Figure 5 Stiffness Distribution 


The strain energy of bending and shear stiffness will be con- 
sidered separately. The bending stiffness is considered first. 
Substituting Equations (3) and (4) into Equation (1) gives the 
strain energy of the bending stiffness represented by one line 
segment i in bay k as 



(5) 


where C = EI^ + |x - x^J (^^q ” ^^p)/(^q ~ ^p) * subscripts 

p and q have been introduced to handle the possibility of a line 
segment extending past the bay limits. Thus x^ is the greater of 


^1 ^ ^q lesser of X 2 or ^+ 1 ' Similarly, EI^ is 

either EI^^ or EIj^ and El^ is either EI 2 or The integration 


is continued for the line segment in adjacent bays, if necessary, 
until the entire line segment has been used. Performing tlie in- 
tegration of Equation (5) yields the strain energy of bending 
stiffness as 
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where 


and 




k.k k,k+l 

^k+l.k+1 



’ V, " 


k 


M, 


k 


f _ 1_ /u2 _ u2\ _ tts. + 

^k.k 2b \“q “pj b2 b3 EI^ 


L El 

f _ —2. _ ^ ^ 3. 

k,k+l b El 

P 


1 El 

f ~ ~ Ati/ ^ 

'k+l,k+l ~ b El 

P 


/El H - El H \/l 

V p q q p// F 

)A 


( 


b » El - El \/L 
q p// p 


H =« X - X, 
p p k 


H = X - »k 


q q 


L = X - X . 

p q p 


( 6 ) 

(6a) 

(6b) 

(6c) 

(6d) 

(6e) 

(6f) 

(6g) 

(6h) 


For constant bending stiffness, l.e., El » £1 , Equations (6a) 

P q 

through (6c) are of Indefinite form. For this case. Integration 
of Equation (5) yields 


f = 1 


l/sEl 

(61) 

k,k 1 

1 q pi 

'/ P 

^k,k+l “ ( 

- h2] 

^ q Pi 


(6j) 

. . , ... » L /El 
k+l,k+l p/ 


(6k) 


where H , U , and L are given above. 
P q P 
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The strain energy of shear stiffness Is considered next. The 
shear stiffness distribution Is similar to that shown in Figure 5 
and is given by the equation 

KAG(x) « KAGi + (x - xi) (KAG 2 - KAGi)/(x 2 - xi). 

Substituting Equations (2) and (4) into Equation (1) gives the 
strain energy of the shear stiffness represented by one line seg- 
ment 1 in bay k as 




where C ■ KAG^ + |x - x^j " ^^p)/ (*q ” ’^p) ‘ subscripts 

p and q have similar meaning as has been previously discussed for 
bending stiffness distribution. Performing the integration of 
Equation (8) yields the same equation as was obtained previously 
for bending stiffness and is repeated here as 


<1 [i;!) 


k,k+l 


Now however* 


As before 


^k+l,k+lj I \ 


^ KAG 
^k.k “ b KAG 


•■k.k+l 


'k+l,k+l - 0 


(KAG^ - ’<AGp|/Lp 


X - X . 

q p 
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For constant shear stiffness, i.e., KAG = KAG , Equation (9a) is 

P <1 

of indefinite form. For this case, integration of Equation (5) 
yields 



L /KAG . 
P/ P 


(9f) 


Equations (6) are evaluated for every line segment of bending 
stiffness and Equations (9) are evaluated for every line segment 
of shear stiffness. Then, application of Castigliano's Theorem 
gives the lateral translation and rotation of panel point k tela-' 
tive to panel point k+1 as 




;ju, 

k 


!l!k 

3Mk 


r* ” 

h,k *k,k+l 


» m 
\ 





( 10 ) 


,'olvlng Equation (10) for and in terms of and and 

substituting into Equation (6) [or (9)] gives the strain energy 
in bay k as 




*^66 68 
(sym) Kqq 


A6, 


A8. 


where 


and 


K 


66 


k+l.k-t-1 

D 


K 


60 


k-fl.k-t-1 


K 


k.k 


60 


^ “ ^k,k ^k+l,k+l " ^Lk+1. 


( 11 ) 


(lla) 

(llb) 

(llc) 


(lid) 
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The restraint at panel point k + 1 is removed by application of 
the transformation 


1-10 


0 0 1-1 


where 

Li ® X, , , - X, . 

k k+1 k 

Using this transformation in Equation (11) gives the final strain 
energy expression as 


k,k k,k+l k,k+n k,k+n+l 




(sym) 


k+1 , k+1 k+1 , k+n k+1 , k+n f 1 

*k+n , k+n ^k+n , k+n+1 


'k+n+1 , k+n+1 I I k+lj 


where 


k.k 6,6 


\,k+l " "‘^6,6 


\,k+n “ ’^60 


^k, k+n+1 “ “^k ^66 " *^60 
\+l,k+l “ *^66 


*k+l,k+n “ "^60 


\+x, k+n+1 “ ^k *^66 ^60 


^k+n,k+n " ^06 


(13h) 
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\+n.k+n+l “ H 

k+n+l,k+n+l tc 66 k 6o cO 

and n is the number of panel points. The kernel matrix in the 
triple matrix product of Equation (13) is the stiffness matrix 
vhich represents the bending and shear stiffness in bay k. 

The stiffness :.atrix for the entire beam is obtained by eval- 
uating Equations (11) and (13) for every bay. Each resulting 4x4 
bay stiffness matrix is added to pi^^vxous 4x4 bay stiffness ma- 
trices at like panel points to form the stiffness matrix (free- 
free) for the entire beam. 

SPECIAL CASE 


For constant bending stiffness (El) and infinite shear stiff- 
ness extending to the bay limits I'x. and x, the stiffness 

loatrix is \ ^ ^ I 


~12 

-12 


( 

1 


12 



(sym) 









where 




\+l 


- X 


k' 


(13i) 

U3j) 
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Subroutine TERPl performs interpolation assuming a linear 
function between known points. The x and y coordinates of the 
known points are given by the elements of {XA} and tlie corre- 
sponding elements in a column of [YA], respectively. Each column of 
[YAj gives the y coordinates of a different set of points. Inter- 
polated y coordinates are calculated at selected x coordinates 
which are given by the elements of IXZ}. These interpolated y 
coordinates are placed in [YZJ. Each column of [Y2] has inter- 
polated values of the respective column of [YA] . Extrapolation 
assuming a linear function is performed when any element of {XZ} 
exceeds the limits of {XA}. 


DERIVATION OF TECHNIQUE 


y 


k 


Given the x,y coordinates of points i and i+1, the coordinate 
at x. is to be found by interpolation assuming a linear function. 



The equation of a straight line is y(x) = ax + b, or in matrix 
notation, y(x) = [x 1] |a| • The coefficients a and b can be 

|b| 

determined because the coordinates x., y. and x..,, V.., are known. 
TU 4 X X x+1 ‘ x+1 

That is. 



a" 

1 

1 

-1 ' 


■yi 

b 

■ ‘‘i+l 

7^i+l 



^1+1. 


from which 
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giving 



The following table gives the correlation between the nomenclature 
of Equation (1) and that used in the Fortran coding in the sub- 
routine. 


Equation (1) 

Fortran Coding* 


XA(1) 

^i+1 

XA(I+1) 

^i 

YA(I) * 

^i+l 

YA(I+1) * 


XZ(K) 


YZ(K) * 


♦Subscript j, denoting different sets 
of points, has been omitted for clarity. 

EXAMPLE 

Consider the following two sets of points oenoted by (1) and 
(2) in the sketch. 
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The coordinates of the points are given by 



-2.' 


-4. 

2. 

iXA} = 

1. 

and r'fAj = 

2. 

-4. 


3. 


1. 

-2. 


(XA) gives the x coordinates of the points in both sets (1) and 
(2)* Column 1 of [YA] gives the y coordinates of the points in 
set (1) and column 2 of [YA] gives the y coordinates of the point?: 
in set (2)* The y values are wanted at x = -1. and x = 7., that 
is, at tXZ) = j-l*j* At X =s -!•, interpolated values of y = 

and y * 0. are calculated from columns 1 and 2 of [YAj , respec- 
tively, using rows 1 and 2 of IXA} and [YAj- At x = 7., extrap- 
olated values of y = -1. and y = 2. are calculated from columns 
1 and 2 of [YA], respectively, using rows 2 and 3 of {XA} and 
[YA]. The final result is , 


[YZ] = 


r-2. 


- 1 . 


0 . 


2 . 


To relate this example problem to a practical problem, con- 
sider IXA) to be the collocation points (panel points) of a vehicle 
and [YA] to be the modal displacements for two modes. Gyros are 
to be placed at stations x = -1. and 7. (i.e., {XZ}). The modal 
displacements ([YZ]) are to be found at these stations* 
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Subroutine TERP2 performs interpolation assuming a diparabolic 
function between known points. A parabolic function is used where 
only three points are available. A diparabolic function in ob- 
tained from the weighted average of two adjacent parabolas and 
will be explained later. The x and y coordinates of the known 
points are given by the elements of {XA} and the corresponding 
elements in a column of [YA], respectively. Each column of [YA] 
gives the y coordinates of a different set of points* Interpo- 
lated y coordinates are calculated at selected x coordinates x^^hich 
are given by the elements of {XZ}. These interpolated y coordi- 
nates are placed in (YZJ. Each column of [YZ] has interpolated 
values of the respective column of [YA]. Extrapolation assuming 
a parabolic function is performed when any element of {XZ} exceeds 
the limits of {XAi. 

DERIVATION OF TECHNIQUE 

The diparabolic interpolation procedure is obtained as follows. 
Because this procedure is dependent upon using parabolas , the 
parabola will be considered first. Given the x,y coordinates of 
points 1» 2, and 3 below^ a parabola is to be fitted to these 
points. 


y 



X 


The equation for a parabola with axis parallel to the y axis is 


or 


where 


or 


y(x) = Ax^ + Bx + C 


y(H) » H 



H 


x - Xj 


X2 - Xi 


H = 


X - X2 


(la) 


X3 - X2 


(lb) 
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is used for ease in later algebraic calculations. The coefficients 
a» b, and c can be determined because the x (or H) and y coordi- 
nates at points 1> 2, and 3 are known* That is» 


y\ 


"1 »i 1 ' 


a 

yi 

= 

R2 H 2 1 


b 



_«3 »3 \ 


c 


from which 


a 


yi 

b 

= M 

Y2 

c 


y3 


where 


H 2 -H 3 


-(H 1 -H 3 ) 


H 1 -H 2 


-(H2-H3HH2+H3) (H1-H3HH1+H3) -(Hi-H2)(Hi+H2) 

_H2H3<H2-H3) -HiH3(Hi-R3) HiH2(Hj-H2> 

(H2-H3) - (H1-H3) + (H1-H2) 


( 2 ) 


therefore 


y(H) = [r2 H 


1] M 


’"I 

V2 . 


ys 


(3) 


For a given set of points, as shown on the following page, a 
parabolic function is used to the left of point 1 and between 
points 1 and 2* Also, a parabolic function is used to the riglit 
of the last point n and between point n-1 and n. Diparabolic 
functions are used between all other points. 



- parabolic - 


■ dlparabolic 


-parabolic - 


Bay 1 (and to the left of point 1) 

I 




The coordinate at is found as follows: 

From Equation (la), 

- Xi)/(x2 - xj) < 

Hi =■ (xi - Xi)/(X2 - xi) = 0 

H 2 = (X2 - Xi)/(X2 - Xi) = 1 

H 3 = (x 3 - Xi)/(X 2 - xi) D. ( 

Using these expressions with Equations (2) and (3), the final 
result Is obtained as 

fl/D 1/(1-D) -1/D(1-D)1 fy] 


y(“k) “ \ “ “k \ 1 |-(1+D)/D -D/(l-D) 1/D(1-D) y;. (4c) 


1 


0 


0 
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The following table gives the correlation between the nomenclature 
of Equations (4a, b, c) and that used in the Fortran coding in the 
subroutine. 


Equations (4a,b,c) 

Fortran Coding* 


XA(m) 


YA(m) * 

«k 

H 


XZ(K) 


YZ(K) * 


*Subscript J, denoting different sets of 
points, has been omitted for clarity. 

Interior Bay i-1 


i+1 





The coordinate yj^ at Xj^ is found as follows. A dlparabollc func- 
tion in bay i-1 of the above sketch is obtained as the weighted 
average of parabolas, A and B. That is, 

y(H) = (1-H) • H • yg (5) 


where 



\ “ *i-l 


as in either Equations (la) or (lb). 


(6a) 
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For parabola A: 

"1-2 ■ (V2 ■ *i-i)/(*i - Vi) 5 c 

*1-1 * (-i-i - - “i-i) ■ » 

“1 ■ (“1 - "i-i)/(«i - *1-1) ■ 1 

From Equations (2) and (3}> 

"^1-2' 

y^(H) - [H^ H 1 ] 

/i _ 

where 

-1/C(1-C) 1/C 1/(1-C) 

- 1/C(1-C) -d+O/C -C/(l-C) 

0 1 0 

For parabola B; 

"1-1 ■ (»1-1 - *i-i)/(*i - Vi) ■ » 
“1 ■ (*1 - ■‘l-l)/(*l - ”1-1) ■ 1 

"1+1 " (Vi ■ *i-i)/(*i ■ *1-1) = “ 

From Equations (2) and (3) 

yg(H) = [h2 h 1] |,^gj 

where 

r l/D 1/(1-D) -1/D(]-D) 

UgJ - |-(l-»-D)/D -D/(l-D) 1/D(1-D) 




(6b) 


(6c) 


1 


0 


0 
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Substituting these expressions for y. and y into Equation (5) 

A D 

and evaluating at point k results in 


where 



l«l<] 


1 

C-D 

D-C 

C(l-C) 

CD 

(1-C) (1-D) 

-2 

2D-C 

1 - 2D + C 

C(l-C) 

CD 

(1-C) (1-D) 

1 

-d+c) 

-C 

C(l-C) 

C 

1-C 


-1 

D(l-D) 

1 

D(l-D) 

0 


0 10 


0 


(6d) 


The following table gives the correlation between the nomenclature 
of Equations (6a, b, c, d) and that used in the Fortran coding 
in the subroutine. 


Equations (6a»b^c,d) 

Fortran Coding* 


XA(I-m) 

i-m 

^l-m 

YA(l-m) * 

H, 

H 

k 


XZ(K) 

1 

^ ^ ■ 1 

YZ(K) * 


*Subscript J, denoting different sets of 
points, has been omitted for clarity. 
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Bay n-1 (and to the right of point n ) 



-L 


The coordinate at is found as follows. 

From Equation (lb) 

\ ■ (\ - Vl)/(*n • -n-l) 

"n-2 ■ (”n-2 ‘ Vl)/(*n ’ \.-l) " 

H 1 “ /x , - X i\//x - X ■ 0 

n-1 \ n-1 ^ 

H ■ /x - X ,\ /ix - X * 1 

n ^ n n-x|/^ n n-l| 

Using these expressions with Equations (2) and (3) » tiie final re- 
sult is obtained as 




-1/C(1-C) 1/C l/(l-c)~ 


V2 

»k ■ K "k 

1/C(1-C) -(l+o/c -C/(l-C) 


^n-1 


0 10 


Jn _ 


The following table gives the correlation between the nomen- 
clature of Equations (7a, b, c) and that used in the Fortran coding 
in the subroutine. 
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Equations (7a,b,c) 

Fortran Coding* 

X 

XA(NXA-m) 

n-m 


1 

3 

YA(NXA-m) * 

»k 

H 


XZ(k) 


YZ(K) ^ 


’^Subscript j, denoting different sets of 
points, has been omitted for clarity. 


example 


Consider the following two sets of points denof-d by (1) and 

( 2 ). 



The coordinates of the points are given by 




• 

- 

10. 


1. 

-]. 

20. 


2. 

-2. 

50. 

and [YA] = 

3. 

-f> . 

70. 


9. 

-2. 

110. 


5. 

•* 

2. 


{XA} gives the x coordinates of the points in both sets (1) and 
(2). Column 3 of [YA] gives the y coordinates of the points in 
set (1) and column (2) of [YA] gives the y coordinates of the 
points in set (2). The y values are wanted at x = S., x = 100., 
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and X = 65., that is, at 


(XZ) 


5. 

100 . 

65. 


At X = 5., extrapolated values of y = 0.375 and y = -0.5625 are 
calculated from columns 1 and 2 of [YAj, respectively, using rows 
1, 2, and 3 of {XA} and [YA]. At x = iOO., interpolated values 
of y = 8- and y = 1.3 are calculated from columns 1 and 2 of [YA], 
respectively, using rows 3, 4, and 5 of {XA} and [YA]. At x = 65., 
interpolated values of y = 7.775 and y = -3.03125 are calculated 
from columns 1 and 2 of [YA], respectively, using rows 2, 3, 4, 
and 5 of {XA} and [YA]. The result is 


[YZj 


0.375 -0.5625 

8. 1.5 

7.775 -3.03125^ ^ 


To relate this example problem to a practical problem, con- 
sider {XA} to be the collocation points (panel points) of a vehicle 
and [YA) to be the modal displacements for two modes. Gyros are 
to be placed at stations x = 5., 100., and 65. (i.e., {XZ}). The 
modal displace lents ([YZ]) are to be found at these stations. 


REFERENCE 

Griffin, J.A.: **A Diparabolic Method of Four-Point Interpola- 

tion.** Journal of the Aerona^itical Sciences^ Vol. 28, No. 2, 
Readers* Forum, February 1961. 
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Subroutine TRAE2 calculates time response additional equations. 

That la 

{Z(t)} - [A] {X(t)} + [B] {X(t)} + {C] lX(t)} + ID] lF(t)} + (E) (1) 

• • • 

vhftr« (X(c}}p {X(t)}» {X(t)}» and (F(c)} are the time response 
and force previously calculated in a time response subroutine, for 
example TRSPl. The TRSPI tape is read in the subroutine TRAE2 to 
calculate {Z(t)}« The coefficient matrices [A], [B], [C]» [0], 
and {£} are supplied or omitted, on option, to this subroutine to 
calculate (Z(t)} which may be shear, bending moment, displacement, 
etc, depending on the choice of coefficient matrices* 



TKANS 


Subroutine TRANS calculates the transpose (inte*^change of rows 


and columns) of a matrix* 
posed, then the result is 

IZI 


If matrix to be trans- 


NCAxNRA 


= lA]^ 


where 


/i = 1, NRA\ 

""ji " (j = 1, NCa) 

NRA is the number of rows of [A] , and NCA is the number of columns 

of [A] • 

Theorem ; The transpose of the product of a set of matrices is 

equal to the product of the transposed matrices taken in 
reverse order. That is, 

(lAltfiJlCl)'^ = [C]"^ [B]'^ [Al^. 

Theorem: The transpose of the sum of a set of matrices is equal 

to the sum of the transposed matrices. That is, 

([A] + [B] + [CJ)’^ = [A]'^ + tB]"^ + [C]^. 

EXAMPLE 


If Che input niacrlx is 


[A] 


2x3 


1. 2. 3. 

4. 5. 6. 


Chen 


^^^3x2 " ^^^2x3 


1. 4. 

2. 5. 

3. 6. 
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Subroutine TRSPl solves the second order differential equation 

[A] (X(t)} + iB] iX(t)i + [Cl iX(t);^ = [D] lF(t) tl) 

for 

lX(t)}» (X(t)}, and (X(t)i 

using a fourth order Runge-Rutta numerical integration procedure 
(Ref 1). Matrices [A], [B], [C], and [D] are input directly to 
this subroutine. For a structural problem, [Aj is tbe mass 
matrix, [B] Is the damping matrix, [C] is the stiffness matrix, 
and [D] is tiie transpose of the vibration mode shapes if the 
equations are a modal representation of the structure. The force 
(F(t)} is calculated internally in the subroutine using linear 
interpolation with [TABT] and [TABF] which are both input to the 
subroutine. As an illustration of the use of [TABT] and (TABFj 
consider the following example. 


2.0 


1.0 

.5 

0 5 10 15 20 

The table (A table is defined in this report as a matrix that may 
have incomplete column data in some rows.) giving the independent 
variable coordinate t is 




[TABT] = 0. 20. 


5. 10. 15. . 


The table giving the corresponding coordinates of the dependent 
variable F Is 
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The following values for lF(t)J will be obtained by linear inter- 
polation at t = 7. 


lF(t = 7)j 


.7 


The Runge-Kutta method of numerical integration has several 
desirable features that may be summarized as follows. The method 
is generally considered to have good convergence qualities. 
Forward integration and iteration procedures can sometimes be 
unstable so that a calculated solution oscillates with rapidly 
increasing amplitude about the true solution. The Runge-Kutta 
method does not seem to be so susceptible to this difficulty. 


DESCRIPTION OF TECHNIQUE 


Matrix [A] must be nonsingular because the highest derivative 
will be calculated from modification of Equation (1). 


U(t)} - [AJ-I [D1 (F(t)} - [A]-l [B] lX(t)) - [A]-l (t) ;X(t); 

Runge-Kutta ^numerical integration is then used to itegrate IX(t)) 
to obtain {X(t)} and then again used to integrate (X(t); to obtain 
{X(t)>. 

REFERENCES 

1* S. Gill: "A Process for the Step-by-Step Integration of Dif- 

ferential Equations in an Automatic Digital Computing Machine," 
Proceedings Cambridge Philosophical Society 47:96-108 (1931). 


(la) 
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Subroutine TRSPIA solves the second order differential equation 

[A] lX(t)} + [B] iX(t)l + [C] lX(t)l - ID] (1) 


for 

{X(t)l, lX(t)i, and lX(t)} 

using a fourth order Runge-Kutta numerical integration procedure 
(Ref 1). Matrices [A], [B], [C], and [D] are input directly to 
this subroutine. For a structural problem, [A] is the mass 
matrix, [B] is the damping matrix, [C] is the stiffness matrix, 
and [D] is the transpose of the vibration mode shapes if the 
equations are a modal representation of the structure. The force 
(F(t)} is calculated internally in the subroutine as follows. 

iF(t)} is size (NF^l) and is obtained considering a (1 - cos)/2 
variation in the coordinate force. Consider the following ex- 
amples. 

EXAMPLE 1: 


Missile Penetrating a Stationary Gust at the Rate - VEL - 



Ifote that the ratio GL/VEL is the period 
of the gust. 
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If FMAG^ is defined to be the maximum normal force at station 
(1). then the 1^^ normal force F(t) varies with time as - 

/ ‘^l\ 

F^(t) - FMAG^ Oa) ll - cos I for 0 - - GL 

(df 1 0. 

F (t) - 0 for / ^ 

(GL^d^ 

lAere, d^ « d^^ - | PP^ - PP^j and 
dj^ - VEL>f(t - START!). 

(assumes PP^ is at beginning of gust at t STARTT*) 


EXAMPLE 2 ; 

In the event a sudden envelopment forcing function (zero lag 
time for all stations) is desired, two alternatives exist: 

1) Supply FMAG and PP as scalar quantities and PP^ = PP, 
This sets the lag time to zero, | PP^ - PP^ * 0| ; or 

2) Supply FMAG as a vector and PP as a vector with all 

stations the same and equal to This sets the 

lag time at all stations to zero. 

The Runge**Kutta method of numerical integration has several 
desirable features that may be summarized as follows. The method 
is generally considered to have good convergence qualities* For- 
ward Integration and iteration procedures can sometimes be unstable 
so that a calculated solution oscillates with rapidly increasing 
amplitude about the true solution* The Runge-Kutta method does 
not aaem to Se so susceptible to this difficulty. 
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DESCRIPTION OF TECHNIQUE 

Matrix [A] must ba nonsingular because the highest derivative 
will be calculated from modification of Equation (1). 

U(t)} - lA)-^ [D] {F(t)} - [A]-l (B] iX(t)} - [A]-l tC] {X(t)}. (la) 

Runge-Kutta numerical integration is then used to Integrate 
{X(t)} to obtain lX(t)} and then again used to Integrate (X(t)} 
to obtain (X(t)}. 

REFERENCES 


1. S. Gill: **A Process for the Step-by-Step Integration of Dif- 

ferential Equations in an Automatic Digital Computing Macl»ine. 
Proceedings Cambridge Philosophical Society 47.96-108 (1951). 
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Subroutine TRSPIB Is a modification of Subroutine TRSPl to use one 
less matrix space so that larger size matrices can be used. The 
second order differential equation to be solved is 

{X(t)l + tB] lX(t)} + tC] (X(t)} = IDJ lF(t)}. (1) 

Comparing this equation with Equation (1) of Subroutine TRSPl shows 
that 


TRSPIB 



[B] 

is 

[A]-l 

[B], 

ic) 

is 

lA]-^ 

IC], 

[D] 

is 

[A]-l 

iDl. 


All other comments in the explanation of Subroutine TRSPl apply to 
this subroutine. 



TRSPIC 


Subroutine TRSPIC is a modification of Subroutine TRSPIA to use 
one less matrix space so that larger s'ize matrices can be used. 

The second order differential equation to be suJved is 

lX(t)} + U] U(t)} [C] {X(t)} = ID] lF(t)J. (1) 

Comparing this equation with Equation (1) of subroutine TRSPIA 
shows that 


TRSPIC 



IB] 

is 

lA]-i 

IB], 

IC] 

is 

IA]-1 

IC], 

ID] 

is 

iA]-l 

[D]. 


All other comments in the explanation of subroutine TRSPIA apply to 
this subroutine. 
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Subroutine TRSP2 sulves the second order differential equation 

[A] U(t)} + [li] tX(t)) + [CJ U(OJ - lU) lK(t)l (1) 

for (X(t)}, (X(t)}, and lX(t)) using a third order Newmark’-Chan- 
Beta numerical integration procedure (Ref I). Matrices [A], iB}, 

[C]» and [D] are input directly to this subroutine* For a struc- 
tural problem^ [A] is the mass matrix, [B] is the damping matrix, 

[C] is the stiffness matrix, and [D] is the transpose of the vi- 
bration mode shapes if the equations are a modal representation 
of the structure. The force lF(t)} is calculated internally in 
the subroutine using linear interpolation with [TABT] and [TABF] 
which are both input to the subroutine. As an illustration of 
the use of [TABT] and [TABFj consider the following example. 


2.0 


0 5 10 15 20 



1.0 

.5 

0 5 10 15 20 

The table (A table Is defined in this report as a matrix that may 
have incomplete column data in some rows.) giving the independent 
variable coordinate t is 

(TABT] • 0. 20. 

5. 10. 15. . 

The cable giving the corresponding coordinates of the dependent 
variable F Is 

[TABF] “2. 2. 

.5 1. .5 . 
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(X}i - 


The following values for iF(l)) will be obtained by linear in- 
terpolation at t = 7. 


(F(t = 7)) 


2 . 


.7J. 


DESCRIPTION OF TECHNIQUE 


Because the difference equation contains disp] .cement terms 
for three consecutive time intervals » it is required to express 
a displacement quantity in terms of initial velocity and initial 
displacement in order to start the procedure. Because [A] must 
be nonsingular in this method, it was decided to modify Equation 
(1) into the form 

{X(t)} + [A)"l [B] lX(t)} + [A]-> tc] {X(t)} = tA]"‘ (D] lF(t)i. 

Thus, only three matrices are needed compared to four matrices 
in Equation (1). The startinp equation is thus given as: 

[S“^ P] iX()} + [S“^ Qj IXo) + ish* ts"‘ D] tKli + (S"* R A"' D] IFio 

where 


[P] 

= [IJ + 

h 

2 

lA- 

B] 

= V'2 - 

• ij) h^ [A 


- (k - 

B) 

h^ 

[A-1 

B] [A- 

•1 C] 

IQ] 

“ h [I] 

- 


- B) 

h3 (A- 

B] tA-^ 

[R] 

- (JS - 

6) 

h^ 

[I] • 

t- (^ - 

3) h^ tA~ 

[S] 

“ m + 

h 

2 

lA- 

B] 

+ lih^ 

[A~l C]. 


[I] is a diagonal of ones (unity matrix) 
h is the integration step size (At). 


After obtaining {X)i, the following difference equation is used 
to obtain the rest of the values of {X}. 
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{X}i^l • rs-l T] {x}^ - [S"l U] 1X}^_J^ 

+ Bh2 lS-‘ A-‘ D1 + (i - 2) IF)^ + 

where 

[T] « 2 [I] - (1 - 2S) h2 [A"i C]. 

[U] « [I] - I [A-l B] + Sh2 [A-1 C]. 

The acceleration and velocity are then obtained using the following 
difference equations. 

■ (i\+ii ■ '“’i ■ ‘^’1) A'** 

‘"Wl = l^ll ♦ I I 

Beta (B) Is the parameter of generalized acceleration and can 
be anywhere between 0 and 

REFERENCES 

1. Chan, S. P., Cox, U. L, and Benfield, W. A.: "Transient 

Analysis of Forced Vibrations of Complex Structural-Mechanical 
Systems." JouiPnal of the Royal Aeronautical Society, July 
1962. 
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Subroutine TRSP2A solves the second order differential equation 

[A] lX(t)} + [BJ {X(t)} -f [C] IX(L)1 = ID) iK(l)l (1) 

for (X(t)}, (X(t)}, and lX(t) 1 using a Lhifd older Neuiwar^-Chau- 
Beta numerical integration procedure (Rei i). Mitiices [A], [B], 

[C], and [D] are input directly to this subroutine. For a stru - 
tural problem, [A] is tlie mass matrix, [B] is the damping matrix, 

[C] is the stiffness matrix, and [D] ‘.3 the transpose of the vi- 
bration mode shapes if the equations are a modal representation 
of the structure. The force lF(t)] is calculated internally in 
the subroutine as follows. 

lF(t)} is size (NF^l) and is obtained considering a (1 - cos)/2 
variation in the coordinate force. Consider tiie following f k- 
amples. 

EXAMPLE 1 : 

Missile Penetrating a Stationary Gust at the Rate - VEL - 



from VEL) 


Wott that the ratio GL/VEL la the period 
of the guat. 
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If FHAG^ Is dsfined Co be Che maximum normal force ac scatlon 
(i), Chen Che 1^^ normal force F(c) varies with cime as - 


F^<C) - FMAG 


F^(c) - 0 


1 <■»> (i - -eir) 


for 0 dj^ < GL 


for 


(diio. 

) GL d^ 


«here, d^ - d^ - ( PP^ - PPj 


and 


d^ - VEL»(C - STARTT). 


Assumes PPi Is at beginning of gust at t = STARTT. 


EXAMPLE 2 ; 

In Che evenc a sudden envelopmenc forcing function (zero lag 
cime for all scacions) is desired, cwo alternatives exist: 

1) Supply FMAG and PP as scalar quantities and PP^ • PP. 
This secs the lag cime to zero, ( PP. - PP « O) ; or 

2) Supply FMAG as a vector and PP as a vector with all 
scacions the same and equal to PP^. This secs Che 

lag time at all scacions to zero. 
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DESCRIPTION OF TECHNIQUE 

Because Che difference equation contains displacement terms 
for three consecutive time Intervals, it Is required tc express 
a displacement quantity In terms of initial velocity and initial 
displacement in order to start the procedure. Because [A] must 
be nonsingular in this method, it was decided to modify Equation 
(1) into the form 

{X(t)} + lAl*l IBJ lX(t)} [Aj-l [C] (X(t)} = lAj'J [D] lF(t)}. 

Thus, only three matrices are needed compared to four matrices 
in Equation (1). The starting equation is thus given as: 

iXJi • IS"^ P] {Xo> + ls~^ Q] (Xq) + ch^- ls~^ A“^ D] iF)i + [S’* R A“^ D] (Flo 

where 


IP] 

» [I] + 1 

[A- 

1 B] 

= (>s - 

i 

r lA-‘ C] 


- (h- &) 

h^ 

lA-i 

B] (A- 

1 C] 


IQI 

« h II] - 

(M 

- tf) 

h^ [A" 

‘ B] 

tA-> BJ 

IRJ 

- (^a - s) 

h2 

[1] 

+ ik- 

d) h 

^ lA-‘ B] 

iS] 

- [I] + 1 

lA- 

1 B] 

+ 3h^ 

lA-1 

c] 

UJ 

la a diagonal 

. of 

ones (unity 

matrix) 


h jS the Integration step size (^t) 

After obtaining {X}), the following difference equation is used 
to obtain the rest of the values of {X}. 

* e>>2 IS-* A- D) ♦ (i - 2) (F)^ * (Fj.jl) 


IT] - 2 [I] - (1 - 2e) h2 [A“i CJ. 
[U] - [I] - I [A-^ B] + eh2 [A-^ C], 


where 
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The acceleration and velocity are then obtained using the following 
difference equations* 

‘^<1+1 • ^ I «*i ^ I «*i+i 

Beta (B) is the parameter of generalized acceleration and can 
be anywhere between 0 and h* 

REFERENCES 

1. Chan, S. P-, Cox, H. L., and Benfield, W. A.: "Transient 

Analysis of Forced Vibrations of Complex Structural-Mechanical 
Systems." Joi4rnaL of the Royal Aeronautical Society^ July 
1962. 
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Subroutine TR3 solves the second order differential equation 

t A ] {X(t)i + t B J (X(t)j + f C J lX(t)} = [Dl lF(t)) (1) 

for {X(t)}» {X(t)}» and (X(t)} in closed form. The results are 
used to solve the matrix equation 

{Z(t)} - [AAl (X(t)} + IBB] lX(t)} 

+ tCC] {X(t)} + [DD] (F(t)} + {EEj (2) 

that may be used to determine dynamic loads or displacements. The 
coefficient matrices [AA], etc are supplied or omitted on option. 

For a structural problem^ t A ] is the generalized mass» t B J 
is the generalized damping, t C J is the generalized stiffness, 
and [D] Is the transpose of the vibration mode shapes. Matrices 
tAj, [Bj, and t C ] are of uncoupled form. The force (F(t)} 
is calculated internally in the subroutine using linear interpola- 
tion with [TABT] and [TABF] which are both input to the subroutine. 

As an illustration of the use of [TABTj and [TABF] consider the 
following example. 



The table (A table is defined in this report as a matrix which 
may have incomplete column data in some rows.) giving the inde- 
pendent variable coordinate t is 


[TABT] - To. 20. 

.5. 10. 
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Th« t«ble giving the corresponding coordinates of the dependent 
variable F Is 


[TABF] 




The following values for {F(t}} will be obtained by linear inter- 
polation at r ■ 7. 


{F(t 


7)} 


r 


2 . 


.7J. 


The matrixes [AA] , [BB], [CC], [DD], and lEE) are input, on 
option, to calculate {Z(t)} which may be shear, bending moment, 
displacement, etc, depending on the choice of coefficient matrices. 

DESCRIPTION OF TECHNIQUE 

Equation (1) may be considered as a set of single degree of 
freedom equations because the mass, damping, and stiffness matrices 
are of uncoupled form. Thus the equation to consider is 

mx(t) + cx(t) + kx(t) = f(t) (3) 


where 


f(t) is a row of [D] (F(t)}. 

Equation (3) is solved by using Laplace transformations. 


In general, I [g(c)J 

£ tx(t)] 
£ [x(t)l 
£ [x(t) ] 


OD 


G(s) = 1 e g(t) dt . 

(4) 

sr 

0 


X(s) 

(4a) 

sX(e) - x(0) 

(4b) 

s^X(s) - sx(0) - x(0) 

(4c) 

F(s). 

(4d) 
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Using the above Laplace transforms vith Equation (3) gives 


X(s) 


F(s) . (s 2a) x(Q) x(0) 

m r(s + b^] (s ^ a)^ + (s + a)- + 


(5) 


where 


a 



(5a) 


u 



<5b) 


The time response x(t) is then found by taking the inverse Laplace 
transforioation of Equation (5) . This can be done conveniently by 
solving each term of the equation separately. Because the force 
may be nonzero at the start time* and because the force is as- 
sumed to vary linearly, the forcing function is broken down into 
a step force and a ramp force. Thus, the response is obtained 
as the sum of four separate solutions. That is. 


x(t) = x(t) due to initial displacement 
+ x(t) due to initial velocity 
+ x(t) due to step force 
+ x(t) due to ramp force. 

In the solution of Equation (5), we will also define 


V 


tij = ^ ^ 


Terms such as sin (bt will use a trigonometic expansion to 

avoiu inaccuracy of combining ip with large values of bt. For 
example. 


if 


ip ■ tan 


rl b 


then 

sin (bt + i|/) ■ ” (a sin bt + b cos bt). 

Ths reepense due to initial displacement, initial velocity, 
s step force, end a ramp force are summarized below. These re- 
sults sre derived in Reference 1. 
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1 * Initial Displacement, x(Q) 

x(t) = — g — e (a sxn bt + b cos bt) 

*x(0) 2 ^ 

x(t) » — ^ ^ e sin bt 

D 

x(0) 2 *at / . , , .V 

x(t) « e (a sin bt - b cos bt) 

D 

For the special case c = k = 0 (i.e., a = b = 0) , the above 
equations give no solution. For this case, 

x(t) = x(0) 

x(t) •= 0 

x(t) =» 0 

2. Initial Velocity. x(0) 

' / ^ x(0) -at , . 

x(t) “ ^ * e sin bt 

D 

x(t) = (-a sin bt + b cos bt) 

b 

x(t) = e l(a^ - b- ) sin bt - 2ab cos bt] 

For the special case c = k = 0, (i.e., a = b = 0) , the above 
equations give no solution. For this case, 

x(t) = x(0) t 

x(t) = x(0) 

x(t) = 0 
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Fs -at^ r 

x/t ^ = — r a * I -a sin bt + b cos bt 
\ s/ mb I s s 

For the special case c = k = 0 (i.e., a *= b = 0) , the above 
equations give no solution. For this case, 
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For the s;>f.cial case c ■ k = 0 (i.e., a = b = 0), the above 
equations give ao solution. For this case, 


■ I; '■r 
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REFERENCES 

1. Wohlen, R. L., and Benfield, W. A.: Closed Foi'm Solution of 

the Second Ordei' Differential Equation, Martin Marietta Corpo- 
ration, Denver Division, Dynamics Section, Memo 108, January 
1967. 
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Subroutine UMAMl calculates a matrix that relates Inertia 
plus applied loads to applied loads for an Inertlally restrained 
(free-free) system. 

A physical Interpretation of this transformation matrix may 
be obtained from the following beam with a single degree of free- 
dom (translation) at each panel point. 



The Inertia forces are defined by 

jFjl - - [A] { 2 } 

where [A] Is the mass matrix oi the beam and (z) is the panel 
point acceleration. 

The panel point accelerations are given by 



( 1 ) 


( 2 ) 


where denotes a reference station that may have any value and 
may or may not coincide with a panel point station. 

Let us define 

[RBM] - {1} 

that Is, the rigid body modes for the beam. The rigid body trans- 
lation mode is given by (1} and the rigid body rotation mode is 
given by - x| . Substituting Equations (2) and (3) into (1) 

gives r .. T 



- - [A] [RBM] 


0 


(4) 
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Fox: loads equilibrium, the sum of the forces and the moments 
about a reference station must be zero. That Is 

[paI + Piij’M 


where |F^| Is the applied external force on the beam. 




i'si - “I 


From Equation (5), noting that 

{1}^ 


|x^ - x| 




- [RBM]^ |Fj.{ = |F^{ 


(5) 


( 6 ) 


Substituting Equation (A) Into (6) 


[RBM] [A] [RBM] 


■R 

LQ J 


[RBM]^ |F^J 


from which 


r** 1 



sf 

R 


.0 ^ 



"] 


-1 


[RBM]* [A] [RBM] [RBM] 


Substituting Equation (7) Into (A) gives the Inertia loads In 
terms of the applied loads as 

{Fll = - tA] [MESS] |F^| 


(7) 


( 8 ) 


where we have defined 

[MESS] = [RBM] [[RBM]^ [A] [RBM]J ^ [RBM]^. 


( 9 ) 



UMAMl— 3"' 


Now the applied plus Inertia Irads can be expressed as 

Pa[ + l"i{ = I"a1 - \h\ 

= tl] - [A] [MESS] |F^| (10) 

[[I] - [A] [MESS]] Is the matrix relating the applied plus inertia 

loads to tb^ applied loads output from the subroutine. It can be 
shown that [[^] ^ [A] [MESS]J is independent of the reference 

station 

Any number of rigid body modes may be used; however, this sub- 
routine is limited to six. 

Several important features of |^[1] - [A] [MESS]] exist. The 
first deals with the triple matrix product 

[[I] - [A] [MESS]]^ [eJ [[I] - [A] [MESS]] 

where ^ J structural influence coefficient matrix of a 

structure restrained in a statically determinant fashion. The 
result of the operation is a free-free structural influence co- 
efficient matrix where the grounding restraints of the structure 
have been removed. It is somewhat helpful to think that the 
grounding restraints of the structure have been replaced by in- 
ertial restraints. A second important feature is that the product 

[^ f ] [h] “ 

where is the free-free structural influence coefficient and 

is the free-free structural stiffness matrix. 

An interesting property of [[I] - [A] [MESS]] is that it is an 

iderrrpotent matrix, that is, a matrix whose produce with itself is 
equal to itself. This can be easily shown as follows. 
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[[I] “ [A] [MESSlj [[I] - [A] [MESS]] 

= [I] - [A] [MESS] - [A] [MESS] + [A] [MESS] [A] [MESS]. 

Now 

[A] [MESS] [A] [MESS] 

= [A] [RBM] [[RBM]^ [A] [RBM]j ^ [RBM]^ [A] [RBi-l] 

[[RBM]^ [A] [RBMl] ^ [K£M]^ 

= [A] [RBM] [[RBM]^ [A] [RBM]] ^ [RBM]^ 

= [A] [MESS], 

•• [[1] - [A] [MESS]] [[I] - [A] [MESS]] = [I] - [A] [MESS]. 

This technique was formulated by David J.ang at Chance Vought 
Aircraft in an informal raeraorondum in 1959 and expanded upon by 
myself and Carl Bodley in I960. 



INITY 


Subroutine UNTTY generates a square matrix with i!i:igunal I'lo- 
ment equal to one and all off-diagonal elements equal to /.ere. 
That is. 


z . , = 0. 


In matrix notation^ 




1 . 


(i = j) 
(i * j) 


where N is the size of [Z] (square) . A synonym for the unit> 
matrix is the identity matrix, thus the usual designation as [1]. 


A matrix is unaltered when multiplied by the unity matrix and 
the process is commutative « In matrix notation. 


[AUU = [I][A] = [A). 



UPDATE 


Subroutine ^TDATE transfers data written by Subroutines VTIAPE and 
YVIAFC £ro« on. or more FORMA tapes onto another (new or existing) FORMA 
tape*. It also provides an option to unload and return the read-stapes and 
drives to the system upon completion of updating from them* 

The selection of matrices to be updated from one tape to the other. 

Is accomplished by specifying the matrix location number (see Subroutine 
VXAPE %irlteup) limits. These limits can Include only one matrix, several 
(consecutive) matrices, or all the matrices (i.e*, the whole tape)* 


^Because Subroutine UPDATE both reads and writes on the ”write” tape, and 
because most computer tape drives do not have sufficient tolerance control 
for a mixed mode (l*e*, read/write) operation, the ’Wite" tape should 
Initially be a disk file* After the update, the disk file may be copied 
onto a tape. 
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Subroutine VCROSS calculates the cross product of two vectors 
in three-dimensional space. By definition, the cross product of 

two vectors A and B is a vector Z whose magnitude is the product 

-> 

of the magnitudes of A and B with the sine of the angle between 
their positive directions; that is» 

-y -y -y 

A B = Z 


|z| » jA| jBl sin (A, B) . 

-y 

The direction of Z is perpendicular to the plane determined by 

A and B and so sensed that the configuration of A, B, Z consti- 
tutes a right-hand system (providing the reference axes them- 
selves are a right-hand system) . 

— ► —y 

In terms of the components of A and B we have. 


"l ■ ^2 S ' “3 '“2 


^^2 ■ ^3 "l ■ “l '’3 


*3 *^2 " ®2 


which is the expanded form of the 



If A X B = 0, then either at 

~y -y 

is zero or A and B are parallel. 


determinant 



least one of the vectors (A, B) 
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The cross product can be represented in matrix notation as 



Note that the first matrix is skew-symmetric. 

Theorem ; Cross multiplication is distributive. That is, 

-> -> -*■ -y -> 

(B+C)=AxB+AxC. 

Theorem : Cross multiplication is not commutative and is some- 

times called anti-commutative since 


A X B = -B X A. 

EXAMPLE 

If the input is A = [-3., 5., 2.] and B = [A., 1., -6.], 

then the output will be 



.90596 
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Subroutine VDOT calculates the dot product of two vectors in 
three-dimensional space • The dot product is also called the 
scalar product or inert product. By definition, the dot product 
— ► — ► 

of two vectors A and B is a scalar equal to the product of the 
magnitude of each vector and the cosine of the angle between 
their positive directions; that is, 

A • B = jA] |b 1 cos (A, B). 

In terms of the components of A and B, we have 

A • B a^ b^ + a^ b^ + a^ b^. 

If A • B - 0, then either at least one of the vectors (A, B) is 
zero or A and B are perpendicular. 

The dot product can be represented in matrix notation as 

A • B = {A}^ {B} 

■[*1 “2 ^ 3 ] 



Theorem ; Dot multiplication is distributive. That is, 
A*(B + C)=A*B+A*C 
Theorem : Dot multiplication is commutative. That is, 

A • B = B • A 
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EXAMPLE 

If the input is A = [-3., 5., 2.] and B = [4., 1., -6.], 

then the output will be 



-.42337 
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Subroutine VMl integrates pressure or weight distribution 
along with concentrated forces or weights to obtain shears and 
bending moments at a set of selected points on a beam. The pres- 
sure or weight distribution and the concentrated forces or weights 
are multiplied by an amplification function (ncrmally the accel- 
eration) • 

The x-stations of the selected points are given in IXVEC). 
These x-stations must be in increasing order. 

The distributed pressure, p(x), is assumed to be piecewise 
linear and is represented by straight line segments as shown in 
Figure 1. 



The x-stations of the end points for the line segments giving the 
distributed pressure are independent of the x-stations in (XVEC). 
The line segments representing the distributed pressure may or 
may not be joined but must not overlap. On any interval where the 
distributed pressure is not defined, the pressure is assumed to 
be zero. The distributed pressure is defined in [DIS]. Each row 
of [DIS] represents one nonvertical line segment. The form of 
each row of [DIS] is [xi X 2 pi p^l where X], pi give the first end 
point and X 2 , P2 give the second end point of a line segment. The 
line segments given by [DIS] must be in increasing order of x. 

The distributed weight may also be given by [DIS]. 

The amplification function is also assumed to be represented 
by straight line segments. All statements given above for dis- 
tributed pressure are applicable. The distributed amplification 
function is defined in [AMP]. 
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The concentrated mass items are defined in [CONC]. Each row 
of [CONC] represents one concentrated mass item and contains: 

1) X 9 the station at which the item is attached to the 

a 

beam; 

2) M . the mass of the item; 

c 

3) X , the center of gravity of the item; and 

eg 

4) the moment of inertia of the item about its own 
c 

center of gravity.. 

These four elements are given in th''. form lx M x . The 

L a c eg c J 

attachment station, x . of an item to the beam must be within the 

a 

panel point limits. 

The calculated shears and bending moments at the selected 
points are placed in {ZV} and {2M}, respectively. 
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Subroutine VMTRl calculates a matrix [Z] that can be used to 
obtain internal shears and bending moments of a beam in terms of 
external forces and moments acting on the beam. The internal 
shears and bending moments » as well as the external forces and 
moments, are at specific locations on the beam. These locations 
are specified by the elements of (PP) which is input to the sub- 
routine. 

EXAMPLE 

Consider the beam below with external forces P. and moments 
acting at the specified locations. 

Pi P 2 P 3 

T,(l lid — iid 

-10 ^20 40 


The locations of the external load application points is given by 


{PP} = 


-10. 

20 . 

L 40 .J 


The internal shears (V) and bending moments (M) at the locations 
specified by {PP} can be obtained in terras of the external loads 
by the following equations: 


J 




(1) 

i=i 



j 

j 


"j 'Zi'j " “i) 

i=l 

i=l 

(2) 
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These equations can be expressed in matrix form for the example 
beam as 



where [Z] is the matrix output from this subroutine. 



WRITE 


Subroutine WRITE writes a matrix of real numbers (a Fortran 
term for numbers with a decimal point) on paper. A group of up 
to ten consecutive elements from a row of tlie matrix are printed 
on each line. If all of the elements of a group are zero, printing 
of this line is suppressed. 

Each matrix printed begins on a new page. On each page of 
printout is the page heading given by Subroutine PAGEHD, the name 
of the matrix, and the row size and column size of the matrix. 

This is followed by the matrix data. On any line of matrix data 
the first integer number is the row number of the matrix elements 
on that line. The second integer number is the column number of 
the matrix element in the first data field. The next group of 
real numbers (up to ten) are the values of the matrix elements. 

This group of matrix elements is given in consecutive column order. 



WRITIM 


Subroutine WRITIM writes a matrix of integer numbers on paper. 
A group of up to twenty consecutive elements from a row of the 
matrix are printed on each line. If all of the elements of a 
group are zero, printing of this line is suppressed. 

Each matrix printed begins on a new page. On each page of 
printout is the page heading given by Subroutine PAGEHI), the name 
of the matrix, and the row size and column size of the matrix. 

This is followed by the matrix data. On any line of matrix data 
the first integer number is the row number of the matrix elements 
on that line. The second Integer number is the column number of 
the matrix element in the first data field. The next group of 
integer numbers (up to twenty) are the values of the matrix ele-- 
ments. This group of matrix elements is given in consecutive 
column order. 
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Subroutine WTAPE writes matrix data at the end of existing 
written matrix data on a FORMA tape (disk is preferred, see below). 
Each set of matrix data consists of two logical records. The first 
record contains the matrix heading (tape identification, location 
number, run number, matrix name, number o^ rov;s of matrix, number 
of columns of matrix, date, and the word "dense"). The second 
record consists of the matrix elements. 

A schematic representation of the tape (disk) is given by the 
following sketch. 


Beginning 

of 

tape (disk) 



El 

H2 

E 2 

• ♦ • 

H 

n 1 

L 

n 

EOT 


where 


= Matrix heading of the i^^ written matrix, 

= Matrix elements of the i^^ written matrix, 

EOT ~ End of Tape. Data written by Subroutine WTAPi£ or 

INTAPE that all FORMA tape subroutines recognize as 
being the end of written data. 

Each vertical line is an end of logical record put on by com- 
puter system’s routines. The tape is written in binary form as 
opposed to binary ceded decimal (BCD) form. 

To find the end of written matrix data, a search is made of 
the matrix headings until the EOT is found. For this reason, a 
"new" tape (disk) must be initialized with Subroutine INTAPE so 
that the tape (disk) contains an EOT. A "new" tape (disk) is de- 
fined to be a tape (disk) for which it is desired to start writing 
matrix data at the front of the tape (disk). Thus, a "new" tape 
(disk) could bv> one with obsolete FORMA matrix data on it as well 
as one that has never been written on by the FORMA system. When 
the EOT is found, a bacitspace operation is done over the EOT, and 
then the current matrix heading, current matrix elements, and a 
new EOT is written. 
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A disk is preferred to a tape for the following reason, be- 
cause of the physical separation of the read and write heads on 
most tape drives there may be tape tolerance problems thus back- 
spacing over the EOT is usually not successful. Instead of ending 
up positioned in front of the EOT, the write head is often posi- 
tioned in front of the previous matrix elements |e^ in the above 

sketchj . The cu^^reiii matrix heading will be written over the pre- 
vious matrix elements. This causes problems later when trying to 
read the records written on the tape. To alleviate this problem, 
it is strongly recommended that all FORMA tape subroutines (INT . ti, 
LTAPE, RTAPE, WTAFE, and UPDATE) use an intermediate device such 
as a disk. At the start of a computer run, the existing tape 
should be copied onto the disk by using computer control cards. 
Likewise, at the end of the run, the disk should be copied back 
onto tape by using computer control cards. 



XLORD 


Subroutine XLORD ^.ontalns a sorting procedure used to arrange 
matrix element locations into ascending order. The associated 
matrix elements are arranged correspondingly. This subroutine 
operates only on matrix elements st<'red In core Subroutine 
XLORD scans the elements to be sorted and selects one of two pos- 
sible methods to be used in sorting. The first method (employed 
primarily when the elements are in a random order) is similai to 
SORTAG by R. C. Singleton*. The second method (employed when the 
elements are broken du^u to tv^o previously ordered groups) is a 
merging procedure. 


*R. C. Singleton: un Effiaifirft Algorithm for Eovting jith Minimal 

Research Meroorandt**^, SRI Project 38753^-132, r>eptember 

1968. 



ZERO 


Subroutine ZERO generates a matrix with each element equal to 
zero« That is» 




0 . 


In matrix notation. 


^^^NRxNC 


0 . 0 . 

0 . 0 . 



/i = 1. NR\ 
\ j = 1 , NC / 


jo, . . . 0,J 

where NR is the number of rows of [Z], and NC is the number of 
columns of tZ]. 


This subroutine is useful in setting a matrix array to zero 
before performing subsequent operations such as matrix assembly 
(ASSEM) or revision/addition of one matrix into another (REVADD). 



ZEROLH 


Subroutine ZEROLH sets the elements below the diagonal of a 
square matrix equal to zero. That is, 

=0. (i > j) 


EXAMPLE 


If [A] Is input to Subroutine ZEROLH as 


tA] 


3x3 


1 . 


2 . 


3. 


4. 5. 6.1, 


7. 


8. 9. 


the matrix output from this subroutine will be 


[A] 


3x3 


1 . 


2 . 


3. 


0. 5. 6.1. 


0 . 


0 . 


9. 



ZEROUH 


Subroutine ZEROUH sets the elements above the diagonal of a 
square matrix equal to zero. That is. 


a 


ij 


= 0 . 


(i ■ j) 


EXAMPLE 


If [A] Is input to Subroutine ZEROUH as 

1- 2. 3. 

4, 5. 6. , 

7. 8. 9. 

the matrix output from this subroutine will be 



[A] 


3x3 


4. 5. 0.1. 


7. 


8 . 


9. 



ZZBOMB 


Subroutine ZZBOMB controls a 'computer run after an error mes- 
sage has been encountered in any of the FORMA subroutines. This 
entails 

1) Printing of subroutine name and error number 
where error occurred; 

2) Printing of all of the core memory used in 
octal form; 

3) Program termination. 



